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ABSTRACT 


In  order  to  better  understand  the  collapse  properties  of 
composite  shells,  an  analysis  of  cylindrical  panels  with  large 
cutout  (11%  and  16.7%)  with  simply  supported  and  unsupported 
(free)  vertical  edges  is  conducted  using  two  finite  element 
programs  (STAGSC-1  and  SHELL).  Cutouts,  their  locations,  boundary 
conditions,  and  panel  width  are  varied  to  see  how  each  affected 
the  collapse  load.  In  addition,  the  shape  and  magnitude  of  the 
residual  strain  due  to  placing  the  cutout  in  the  panel  is 
evaluated . 

It  is  found  that  the  panel  collapses  at  lower  load  (a 
reduction  of  50-80%)  when  a  cutout  is  introduced,  if  the  vertical 
supports  are  removed,  or  if  the  panel  width  is  reduced.  The  panel 
collapses  at  higher  loads  (an  increase  of  3-10%)  when  eccentricity 
(17%-25%)  in  cutout  location  is  introduced.  Initial  imperfections 
in  the  solid  panels  create  a  reduced  panel  stiffness  which  reduces 
the  collapse  load.  The  new  version  of  STAGS  (1986)  is  seen  to 
accurately  predict  panel  response  in  terms  of  load  and 
displacement.  SHELL  incorporates  through-the-thickness  shear 
which  reduces  panel  stiffness  and  gives  closer  approximation  of 
experimental  results  than  STAGS.  It  was  determined  that  if  the 
cylindrical  composite  panel  undergoes  rotations  of  greater  than 
15-17**,  or  displacements  larger  than  five  times  the  panel 
thickness,  transverse  shear  effects  must  be  addressed  in  the 
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analytical  modeKs).  This  requires  using  the  SHELL  code  to 
accurately  predict  panel  response  for  axial  compression.  If  the 
panel  does  not  undergoe  such  rigorous  displacements  and/or 
rotations,  then  STAGSC-1  code  will  give  accurate  predictions.  The 
STAGSC-1  code  is  adequate  for  all  the  models  considered,  but  does 
not  incorporate  a  transverse  shear  strain  model.  Hence,  the 
rational  for  incorporating  SHELL  within  this  research  effort. 
Residual  stress  does  exist  when  a  cutout  is  introduced  to  the 
panel .  This  stress  is  significant  ( w  750  psi )  but  does  not 
dominate  the  panel’s  response  when  compared  to  the  effect  of 
geometric  imperfections. 


INVESTIGATION  OF  COLLAPSE 
CHARACTERISTICS  OF  CYLINDRICAL  COMPOSITE 
PANELS  WITH  LARGE  CUTOUTS 


1.  INTRODUCTION 

1.1  Background 


Large,  lightweight  shell  structures  are  used 
extensively  in  aircraft  for  both  strength  and  stability. 
These  shell  structures  are  greatly  influenced  by  either 
reinforced  or  non-reinforced  cutouts,  cutout  eccentricity, 
and  mixed  boundary  conditions.  Due  to  these  influences  and 
with  recent  advances,  a  vast  number  of  these  shells  are  made 
from  composite  materials.  These  composite  materials  are 
desired  primarily  due  to  the  ability  of  the  designer  to 
tailor  the  material  properties  in  a  structure  to  withstand 
the  load  environment.  This  material  tailoring  leads  to,  in 
general,  a  lighter  structure  than  similar  metal  structures 
thereby  reducing  the  overall  weight  of  a  system;  a  critical 
factor  in  aerospace  design.  Therefore,  from  a  practical 
point  of  view,  there  is  a  need  to  further  research  into 
composite  shell  structures  and  the  effect  of  cutouts,  cutout 
eccentricity  and  boundary  conditions  on  the  stability  of 


shell  structures . 


The  basic  area  of  work  on  shells  is  very  broad.  A  good 
portion  of  this  work  deals  with  geometrically  linear  problem 
solving.  Ugral  [1]  shows  a  basic  layout  of  linear  shell 
theory.  Other  works  such  as  Saada  [2] ,  Sanders  [3]  and 
Leissa  [4]  deal  with  shells  from  a  nonlinear  point  of  view. 

Very  few  references  can  be  found  dealing  with  buckling 
of  composite  shell  structures  (such  as  panels  or  plates) 
under  axial  compression  with  cutouts.  In  contrast,  numerous 
studies  have  been  carried  out  on  isotropic  plates  and 
cylindrical  panels  without  cutouts  [5-10] .  These  studies 
center  around  various  boundary  conditions  and  their  affect 
on  the  buckling  load.  If  one  considers  the  area  of  axially 
loaded  isotropic  panels  with  cutouts  [11-15] ,  mainly,  panels 
with  centered  circular  cutouts  were  analyzed  numerically  by 
finite  difference  techniques  to  determine  the  buckling  and 
bifurcation  loads. 

As  mentioned  earlier,  composite  panels  are  an  important 
part  of  aerospace  structures  due  to  their  excellent 
strength-to-  weight  and  stiff ness-to-weight  ratios. 

Composite  cylindrical  panels  without  cutouts  have  been 
studied  [16-20] ,  but  very  little  work  can  be  found  in 
dealing  with  composite  panels  with  cutouts.  Janisse  [21] , 
Lee  [22] ,  Hermsen  [23] ,  Dennis  [24]  ,  and  Knight  and  Starnes 
[30-31]  are  authors  who  have  studied  the  problem  of  cutouts 
in  composite  cylindrical  panels.  Janisse  [21]  found  that 


the  collapse  characteristics  of  composite  cylindrical  panels 
are  dependent  upon  ply  lay  up  and  size  of  the  cutout.  Lee 
[22]  studied  the  effects  of  the  cutout  aspect  ratio  and 
concluded  that  as  the  surface  area  of  the  cutout  increased 
the  buckling  load  decreased  with  the  effects  of  nonlinearity 
becoming  more  pronounced.  Hermsen  [23]  concluded  that 
eccentricity  of  small  cutouts  will  reduce  the  buckling  load 
and  that  the  larger  cutouts  induce  nonlinear  response 
throughout  the  entire  loading  range.  Dennis  [24]  developed 
a  two-dimensional  shell  theory  completely  described  by 
orthogonal  curvilinear  coordinates.  This  theory  encompasses 
large  displacement  and  rotational  geometric  nonlinearity  for 
small  strain  situations.  Additionally,  the  theory  includes 
a  parabolic  transverse  shear  stress  and  strain  distribution 
through  the  shell  thickness.  Linneman  [25]  has  an  excellent 
presentation  of  the  theoretical  development  of  anisotropic 
cylindrical  shell  theory  and  the  corresponding  equations  of 
motion  and  appropriate  boundary  conditions. 

Uith  the  ever  increasing  capability  of  the  modern 
computer,  many  shells  are  analyzed  numerically.  Most  of  the 
recent  computer  work  is  finite  element  based.  Cook  [26]  has 
a  very  well  written  section  on  nonlinear  solution  techniques 
while  the  Structural  Analysis  for  General  Shells  (STAGS) 
theory  manual  [27]  and  Zienkiewicz  [28]  have  sections  that 
deal  specifically  with  shells  and  geometrically  nonlinear 


finite  element  solutions. 

Other  references  are  concerned  with  experimental 
research  on  composite  shell  structures.  This  is  of 
particular  importance  in  order  to  validate  closed  form  and 
finite  element  solutions.  References  in  this  area  include 
works  by  Lee  [29] ,  Tisler  [32-35] ,  and  Egan  [36-37] . 

This  thesis  addresses  all  of  these  aspects  of  composite 
shells  and  includes  three  aspects  that  are  not  well 
documented.  These  three  aspects  incorporate  nonlinear 
analysis  with  experimental  testing.  One  aspect  relates  to 
the  comparison  of  the  shell’s  response  to  axial  loading  when 
the  vertical  edges  are  completely  unsupported  as  compared  to 
simply  supported  edges,  thereby  increasing  large 
geometric  displacement  and  rotational  nonlinearity.  Another 
aspect  is  investigating  the  eccentricity  (17-25%)  of  large 
cutouts  (11-17%  reduction  in  surface  area)  and  the  combining 
this  effect  with  various  boundary  conditions.  The  final 
aspect  relates  to  the  experimental  determination  of 
residual  stress  generated  by  removing  large  cutouts  and 
analytically  determining  the  effect  this  residual  stress  has 
on  the  buckling  characteristics.  There  is  very  little 
documentation  on  laminated  cylindrical  shells  in  the  areas 
considered  by  this  thesis. 


1.2  Objectives 


The  major  purpose  of  this  thesis  is  determining  if  a 
nonlinear  finite  element  analysis  of  a  laminated  cylindrical 
shell  produces  results  that  correlate  well  with  results  from 
experimentally  tested  panels.  The  program(s)  used  to  do 
this  analysis  are  ( 1 )  the  Structural  Analysis  for  General 
Shells  (STA6SC-1)  computer  program  (1986  version)  and  (2) 

Dennis’  [24,40-41]  two-dimensional  theory  (SHELL)  that 
incorporates  a  parabolic  transverse  shear  strain  energy 
distribution  through  the  thickness.  The  shell(s)  are 
axially  loaded  with  an  external  distributed  load  along  the 
top  surface. 

The  second  purpose  ir  investigating  the  effect  of 
various  boundary  conditions,  simple  supported  vs  unsupported 
vertical  edges,  (see  Figs  1.'’  and  1.2),  and  eccentricity  of 
the  cutout  location  to  the  shells’  collapse  characteristics 
(Figs  1.3  and  1.4),  and  determine  what  effect  transverse 
shear  strain  energy  has  during  large  geometric  displacement 
and/or  rotations. 

A  third  purpow  aS  comparing  analytical  and 
experimental  results.  The  final  purpose  is  a  residual 
stress  analysis.  This  analysis  is  concerned  with 

determining  the  magnitude  and  shape  of  the  stress  when  removing  the 
cutout,  then  determining  the  effect  of  this  stress  to  the 


{u=R^=  free 
v=w=R  =R  =0 

V  w 


1  Panel  with  Simply  Supported  Vertical  Edges 
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Fig  1.2  Panel  with  Unsupported  Vertical  Edges 
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Fig  1.4  Panel  with  4“  Cutout,  Offset  1" 
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shell’s  collapse.  The  results  are  used  to  better  understand 
the  effects  of  large  displacements  and  rotations  on  a  finite 
element  analysis. 

1.3  Scope 


A  total  of  tMelve  models  are  run  on  STAGS  and  six 
models  are  run  on  SHELL  to  determine  displacements  due  to 
an  external,  axially-applied  distributed  load.  The  models, 
developed  from  previous  work  [35,41] ,  included  the 
material  properties  of  the  experimental  test  specimens  and 
incorporated  a  modified  Riks  technique  based  on  the 
"constant  arc  length"  method  [39] .  The  models  are  run  with 
a  nonlinear  analysis  using  STAGSC-1  and  SHELL  with  linear 
analysis  runs  for  comparison  purposes. 

Next,  the  results  are  compared  to  experimental 
results.  The  experimental  results  are  obtained  by  the 
Flight  Dynamics  Laboratory,  Wright  Research  &  Development 
Center,  Wright-Patterson  AFB,  Ohio. 

Finally,  experimentally  derived  residual  stress 
parameters  are  incorporated  into  one  nonlinear  analysis 
model  to  compare  residual  stress  effects  on  nonlinear 
collapse  characteristics. 


2.  THEORETICAL  DEVELOPMENT 


The  first  step  in  the  theoretical  development  for  this 
thesis  is  the  derivation  of  a  displacement  field  for  a 
general  shell  based  upon  nonlinear  kinematics.  This  is  done 
to  give  the  reader  an  understanding  of  the  complex  nonlinear 
relationships  that  govern  a  arbitrary .cylindrical ,  composite 
shell’s  response  to  axial  loading.  Next,  the  theory 
governing  STAGSC-1  analysis  is  developed  for  a  plate 
element.  This  development  includes  kinematics,  classical 
laminated  plate  theory,  and  the  solution  techniques.  Next, 
the  theory  governing  SHELL  analysis  is  developed  similarly 
for  a  cylindrical  shell.  The  work  developed  subsequently  is 
a  review  of  many  references  on  shell  theory.  It  is 
presented  to  the  reader  for  completeness  only.  The  writer 
does  not  intend  to  imply  originality. 

2.1  Derivation  of  General  Nonlinear  Shell  Strain 
Displacement  Equations 

Since  this  thesis  is  concerned  with  three  dimensional 
curvilinear  shells,  a  basic  understanding  of  the  underlying 
kinematics  (in  orthogonal  curvilinear  coordinates)  is 
necessary.  Therefore,  a  brief  presentation  is  made  of  the 
derivation  of  the  strain  tensor  giving  the  exact  relations 


in  orthogonal  curvilinear  coordinates.  These  relations  are 
presented  in  several  elasticity  texts,  see  [2]  for  example. 


Consider  a  point  M  of  Fig  2.1  that  is  located  in  3-D 
space  by  the  position  vector,  r.  The  point  M  has  Cartesian 
coordinates,  x^,  as  shown  in  Eqn  (2.1). 

r  =  X.  I.  i  =  1,2,3  (2.1) 

1  i 

where , 

X  =  Cartesian  coordinates 

i 

I  =  Cartesian  basis  vectors 

i 

A  summation  convention  on  repeated  indices  applies 
unless  otherwise  stated,  and  barred  quantities  refer  to 
vectors . 

The  point  M  can  also  be  represented  with  curvilinear 
coordinates  y  .  The  Cartesian  coordinates  are  related  to 

i 

the  curvilinear  coordinates  through  transformations  of  the 
form  shown  in  Eqn  (2.2): 


X  *  X  (y  .y  fy  )  =  x  (y  ) 

i  1  1  Z  3  i  j 

y  =  y(x,x,x)  =  y(x) 

i  i  l  2  3  i  j 


(2.2) 


Basis  vectors,  a^,  of  the  curvilinear  system  are  found 
by  taking  the  differential  of  the  position  vector,  r. 


Fig  2.1 


Point  M  Located  in  3— D  Space  by  Position  Vector 
Vector  dr  has  length  ds. 
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From  Eqn  (2.1) 


df 


dx  i 

i  i 


(2.3) 


The  length,  ds,  of  the  infinitesimal  line  segment  MN  (see 
vector  dr  of  Fig  2.1)  is  then  given  by. 


(ds)^  =  dr  •  dr  (2.4) 

The  length  ds  is  independent  of  coordinate  systems  and 
therefore  from  Eqn  (2.2), 


dr  =  a  dy  +  a  dy  +  a  dy 

1  1  2  2  3  3 


(2.5) 


where 


a  = 


d  r 


r , 


a  y. 


A  curvilinear  basis  vector,  a^,  is  tangent  to  the 
coordinate  line.  From  Eqns  (2.4)  and  (2.5)  we  can  write. 


(ds)"'  =  (a 


a  )dy  dy 
j  i  J 


=  g 


ij 


dy  dy 
i  J 


(2.6) 
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where 


The  elements  of  g  form  a  symmetric  tensor  called  the 

» j 

metric  tensor  that  links  the  two  coordinate  systems  and 
through  the  invariant  property  of  length.  The 
coordinate  system  is  called  orthogonal  if  its  metric  tensor 
is  diagonal,  i.e.,  when  ~  ^  i  J.  This  is  assumed 
from  this  point  on. 

Next,  consider  the  infinitesimal  line  segment,  MN,  of 

length  ds  now  embedded  in  a  differential  volume  element,  see 

Fig  2.2.  This  differential  volume  is  linearly  transformed, 

i.e.,  deformed,  to  a  new  configuration  where  the  line 

segment  is  now  of  length  ds  ,  and  whose  transformed 

coordinate  system  has  a  metric  tensor  G  .  As  a  result  of 

ij 

deformation,  the  line  segment  MN  of  Fig  2.2  moves  to  M  N 
represented  by  the  displacement  vector,  u.  By  subtracting 
the  original  and  deformed  squared  lengths  of  the  line 
segment,  the  Green’s  strain  tensor,  7^^,  is  defined  as  shown 
in  Eqn  (2.7). 


(ds*)^  -  (ds)^  =  27^jdy^dyj  (2.7) 


Previously  defined  quantities  then  give 


=  G.,  -  9..  -  a  ‘U,.  +  a.*u,.  +  u,.*u. 

ij  ij  »J  i  J  J  1  1  j 


(2.8) 


The  physical  strains. 


are  then  found  from. 


e 


ij 


h  h 
»  j 


(2.9) 


In  Eqn  (2.9),  the  h^  are  called  scale  factors  and  are 
2 

defined  by  g  =  h  (no  sum)  and  the  y  are  shown  in  Eqn 

ij  i  'ij 

(B.IO)  where  the  u^  are  the  coordinates  of  the  displacement 
vector,  u.  In  the  general  large  strain  case,  the  e^^(no 
sum)  are  related  to  the  elongations  of  the  fibers  of  the 
differential  volume  element  and  the  (i#j)  are  related  to 
the  shears  (the  difference  from  ninety  degrees  the 
originally  perpendicular  fibers  are  oriented  to  after 
deformation).  For  the  case  of  small  strains  (e  <  .04)  [41], 
the  e^^(no  sum)  are  the  elongations  and  the  e^^(i#j)  are  the 
shears.  That  is  to  say,  for  large  strains,  the  c^^have  no 
physical  meanings;  but  for  small  strains,  the  e^^have  the 
engineering  definitions.  The  elongations  and  shears  are 
identically  zero  for  rigid  body  displacements  and  rotations 
of  the  body  under  loading  [41].  In  this  way,  there  are  no 
theoretical  limitations  on  the  magnitudes  of  displacement 
and  rotation  elements  of  which  the  body  can  undergo  [41] . 
One  important  limitation  of  Eqn  (2.10)  is  that  it  applies 


2-6 


only  to  small  strain  quantities.  Therefore,  if  large 


•  strains  exist,  Eqn  (2.10)  is  no  longer  valid.  This  theory 

applies  only  for  large  displacements  and/or  moderately  large 
rotations  with  small  strains  (e.g.  no  plasticity). 


Fig  2.2 


Segment  MN  Deforms  to  Through  Displacement 

Vector  u. 
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2.2  STA6SC-1  Theory 


The  Structural  Analysis  for  General  Shells  (STAGS) 
computer  program,  developed  by  Lockheed  Palo  Alto  Research 
Laboratory,  was  first  operational  in  1967.  Its  primary 
purpose  is  the  analysis  of  thin  shelled  structures.  From 
1967-1976  the  program  was  based  on  the  finite  difference 
method.  In  1979  a  new  version  of  STAGS,  STAGSC-1 ,  was 
released  and  is  based  entirely  on  the  finite  element  method 
[42] .  The  1986  STAGSC-1  VAX  Computer  version  is  used  in 
this  thesis.  Of  the  many  capabilities  in  the  STAGS  program, 
the  primary  one  used  herein  is  the  nonlinear  static  analysis. 
The  basic  finite  element  used  in  STAGSC-l  is  the  flat  plate 
element  that  facets  the  shell  surface  to  approximate  its 
curvature.  A  look  at  the  theory  for  this  flat  plate  element 
and  the  theory  used  in  STAGS  is  presented  in  the  following 
sections,  based  primarily  upon  Egan’s  work  (see  [37]  ). 

Uhile  this  may  not  be  the  exact  theory  used  in  STAGSC-1,  the 
following  presentation  is  an  attempt  through  researching 
documentation  to  give  some  insight  into  the  internal  make-up 
of  the  STAGS *8  code. 


2.2.1  Nonlinear  Shell  Theory 

A  review  of  Sanders’  nonlinear  shell  theory  is 
presented  so  that  the  reader  may  understand  strain- 
displacement  relations  in  a  shell.  This  presentation  Mill, 
in  the  subsequent  section,  provide  an  insight  to  the  theory 
used  in  STAGSC-1 . 

The  classical  thin  shell  theory  derived  by  Sanders 
assumes  that  the  shell  is  thin,  the  middle  surface  strains 
and  rotations  are  small,  and  the  displacements  away  form 
the  surface  are  restricted  by  the  Kirchhoff-Love  hypothesis 
[3] .  Uith  these  assumptions  in  mind  the  equations  for  the 
mid-surface  and  bending  strains  are  [44] : 
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where  u  and  v  denote  displacements  tangent  to  the 
mid-surface  and  m  denotes  displacement  perpendicular  to  the 
shell’s  mid  surface,  R  and  R  are  the  principal  radii  of 

X  y 

curvature,  ,  <tf  ,  and  0  are  in-plane  rotational  terms,  and 
X  y 

a  and  a  are  metric  tensor  coefficients.  The  rotational 

X  y 

terms  are  given  by  [44] : 


(2.13) 


The  metric  coefficients,  written  in  a  form  similar  to  Saada 
[2] ,  are: 


(2.14) 


where  E  is  the  distance  from  the  mid-surface  and  E  and  G  are 
3 

functions  of  the  shell  geometry  (not  material  constants). 

For  simplicity,  a  thin  cylindrical  shell  is  considered 
for  the  remainder  of  the  development  (see  Fig  2.3).  The 
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reference  surface  for  the  shell  is  the  mid-plane.  For  this 
cylindrical  shell  we  let  R  -»<»,  R  =  R.  E  =  z,  and  assume 

X  y  ’3 

that  z  «  R.  For  this  particular  case  Sanders’  equations 
(Eqns  (2.1)  and  (2.2))  reduce  to  [17]: 
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Now,  Eqns  (2.15)  and>^2.16)  can  be  written  with  the 
Kirchhoff-Love  hypothesis  to  obtain  the  full  strain 
expression: 
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2.2.2  Strain-Displacement  Relations 


As  stated  earlier,  STAGS  uses  flat  plates  to 
approximate  the  curved  surface  of  a  shell.  These  plate 
elements  are  thin,  therefore  a  state  of  plane  stress  can  be 
assumed  with  y  ,  y  ,  e  ,  and  a  equal  to  zero,  and  the 

xz  yx  z  z 

in-plane  displacements,  u  and  v,  as  well  as  the  normal 
displacement,  w,  functionally  depending  on  only  two  space 
variables  [26].  As  in  Sanders*  equations  for  a  thin  shell, 
STAGS  uses  the  Kirchhoff-Love  hypothesis  for  strains  away 
from  the  mid-plane.  The  complete  derivation  of  the  nonlinear 
kinematics  relations  is  given  in  Appendix  A. 

If  the  Kirchhoff-Love  hypothesis  is  considered  for  out 
of  plane  strain  terms  (see  Appendix  A)  and  combined  with  the 
in-plane  strain  terms  from  Eqn  (2.20),  the  resulting 
expressions  for  strains  in  the  plate  are  given  by: 
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(2.18) 


Eqn  (2.18)  shows  the  nonlinear  kinematic  equations 

which  appear  to  be  used  in  STAGSC-1 .  These  kinematic 
relations  allow  for  large  displacements  and  moderate 


rotations  (due  to  Kirchhoff-Love  hypothesis)  and  will  be  used 
in  the  following  section  to  demonstrate  the  general  form  of 
the  tangent  stiffness  matrix. 


2.2.3  Constitutive  Development 

As  was  stated  in  Section  2.2.2,  a  state  of  plane 
stress  is  assumed  for  the  plate  element.  The  following 
contracted  notation  in  Table  2.1  is  used  for  the 
constitutive  development  where  Fig  2.4  demonstrates  the 
relationship  between  the  subscripts. 

Table  2.1  Contracted  Notation  for  STAGS 


stress 

strain 

plate 

explicit {contracted 

explicit {contracted 

coordinates 

o  a 

11  1 

Si 

X  -»  1 

a  a 

22  2 

S2  S 

y  -*  2 

a  a 

33  3 

e  e 

33  3 

z  -»  3 

shear : 

a  a 

23  4 

2e  e 

23  4 

y-z  -»  4 

a  a 

13  5 

2e  c 

13  5 

x-z  5 

a  a 

12  6 

2e  £ 

12  6 

x-y  -»  6 
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The  research  in  this  thesis  is  limited  to  material 
response  in  the  linear  regime,  but  complexity  arises  due  to 
the  directional  response  of  the  oriented  plies  in  the 
laminate.  The  constitutive  relations  between  stress  and 
strain  for  a  laminate  of  arbitrarily  oriented  transversely 
isotropic  plies  is  developed.  These  relations  are  then 
summed  over  the  lamina  thickness  of  the  directional 
constitutive  equation  for  each  ply,  to  arrive  at  the  total 
laminate's  effective  stress-strain  relation.  See  Appendix  D 
for  complete  derivation  based  upon  [45] . 


Fig  2.4  Fiber  Reinforced  Lamina  Definitions 


From  basic  elasticity  relations,  for  isotropic  material 
under  purely  extensional  load,  stress  (a)  and  strain  (c)  are 
related  as: 


=  Ee 


(2.19) 


with  E  being  Young’s  modulus.  For  a  general  anisotropic 
material.  Young’s  modulus  can  differ  with  different  load 
orientations,  so  the  equation  expands  to  [45]: 
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where  is  the  stiffness  matrix  of  36  independent  terms 

defining  the  stress-strain  relationship  for  loading  in  the 

ith  direction.  Uhen  the  area  of  study  is  restricted  to  the 

energy  conserving  elastic  regime,  the  matrix  is  symmetric 

(S  -  S  )  resulting  in  21  independent  terms.  In  the  case 
12  21 

of  fiber-reinforced  composites,  shown  in  Fig  2.4,  the  three 
mutually  orthogonal  planes  of  symmetry  decouple  shear 
strains  from  normal  stresses  and  vice  versa,  and  shear 
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stresses  are  decoupled  from  shear  strains.  This  is  the 
definition  of  an  orthotropic  material,  which  has  only  nine 
independent  terms  [45] : 
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Since  the  material  responds  equally  to  any  direction  of 
load  in  the  plane  perpendicular  to  the  fiber  longitudinal 
axis,  so  the  2  and  3  subscripts  on  S  are  interchangeable. 
This  behavior  is  termed  transverse  isotropy,  and  it  further 
reduces  the  number  of  independent  terms  to  seven.  Jones 
defines  these  stiffness  terms  with  respect  to  the 
engineering  constants  E  (Young’s  Modulus)  and  u  (Poisson’s 
ratio)  [45]  (refer  to  Appendix  E,  Eqn  (E.2)). 

Due  to  the  thinness  of  the  plies  in  the  laminate,  the 
assumption  of  plane  stress  is  made  (a  -  o  -  o  =0)at 

3  4  5 

this  point.  It  should  be  noted  that  the  limitation  to  this 
assumption  is  plane  stress  is  no  longer  valid  for  thick 
plies.  Using  a  simplistic  criterion  that  t  ^  l/2n  inches 
(see  [54] )  as  the  point  when  plane  strain  applies,  one  can 


easily  determine  that  the  laminate  (t  =  .004")  is  much  less 
than  the  criteria  and  the  plane  stress  assumption  is  valid. 
Solving  for  in  Eqn  (2.21)  after  applying  this  assumption 
yields  [45] : 


=  -  (S  /S  )£  -  (S  /S  )e 

13  33  1  23  33  2 


(2.22) 


Applying  the  plane  stress  assumption  to  Eqn  (2.24)  and  using 
the  relation  of  Eqn  (2.21)  to  eliminate  £  produces  the 

3 

lamina  constitutive  relation: 
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where  the  are  the  reduced  stiffness  coefficients  related 
to  the  S^j’s  by  [45]: 
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ij 


ij 


( S  S  )/S 
13  j3  33 


(2.24) 


where  is  defined  in  Appendix  D.  Finally,  in  order  to 

analyze  the  stack  of  plies,  they  must  all  be  referenced  to  a 
global  axis  system  and  their  effects  summed  [45] : 


(2.25) 
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where  [T]  is  defined  in  Appendix  E.  Uith  this 
transformation  the  constitutive  relations  become: 


{•■I- 


(2.26) 


with  the  transformed  reduced  stiffness,  Q.  ,  being  defined 


ij 


in  Appendix  D. 


2.2.4  Derivation  of  the  Tangent  Stiffness  Matrix 


The  tangent  stiffness  matrix  is  a  nonlinear  stiffness 
matrix  used  in  the  modified  Newton-Raphson  method  for 
solving  a  nonlinear  set  of  equations.  There  are  other 
techniques  for  solving  these  equations,  but  STAGSC-1  uses 
the  modified  Newton-Raphson .  Therefore,  the  derivation  of 
the  tangent  stiffness  matrix  is  of  interest  for  this  thesis. 
The  tangent  stiffness  matrix  is  derived  for  a  flat  plate 
element,  a  STAGS  type  element,  in  a  general  way  without 
reference  to  any  specific  shape  functions.  The  specific 
shape  functions,  and  their  derivation,  for  the  STAGS  element 
is  shown  in  [43] .  It  must  also  be  pointed  out  that  STAGS 
uses  an  isoparametric  formulation,  which  requires  the 
Jacobian  for  the  element,  whereas  this  formulation  is  given 


2-21 


in  general  coordinates  to  show  the  reader  the  steps  involved 
in  formulating  a  nonlinear  stiffness  matrix.  A  complete 
derivation  of  the  tangent  stiffness  matrix  is  shown  in 
Appendix  F  and  in  [37] . 

To  start  the  derivation  (see  Appendix  F  for  the 
complete  derivation),  a  form  of  the  tangent  stiffness  matrix 
is  found  from  the  energy  expression.  Consider  (vi)f  the  sum 
of  the  external  and  internal  generalized  forces  (Appendix  C, 
Eqn  (C.14)),  which  is  given  as  [28]: 


where , 

(a)  =  vector  of  stresses 
( f )  =  nodal  forces  of  the  element 
In  an  equilibrium  state,  (ip)  (Eqn  (2.27))  will  equal  zero. 
The  strain  displacement  relations  can  be  written  (to  include 
the  [B]  matrix)  as  [37]: 

(e)  =  [L][N](a)  =  [B](a)  (2.28) 


where , 

[L]  =  differential  operator  matrix 
[N]  -  element  shape  functions 
( a )  -  nodal  displacement  vector 


In  the  case  of  a  nonlinear  stiffness  matrix  (i.e. 
tangent  stiffness  matrix),  [B]  is  redefined  as  [26]: 

[B]  =  [Bo]  +  [Bl]  (2.29) 


where , 

[Bo]  =  constant 

[Bl]  =  functions  of  displacements 

In  order  to  use  the  Newton-Raphson  method,  a  relation 
between  d(a)  and  (see  Section  2.2.7)  must  be  found  [37]. 

Taking  the  variation  of  Eqn  (2.27)  with  respect  to  5(a) 
gives  the  relation  needed  (see  Appendix  F,  Eqn  (F.4)). 
Incorporating  the  elastic  stress  strain  relations  into  the 
variational  equations  and  defining  stiffness  matrices  based 
on  the  appropriate  terms  yields  [27] : 


[Ko] 

[Kl] 


f  [Bo]  ”^[0]  [Bo]  dV 
V 

J  [[Bo]^[D][Bo]  +  [Bi]^[D][Bl] 

+  [Bi][D][Bo]l  dV 


dV 


(2.30) 


[Ka]  6(  a  ) 


where 


[Ko]  =  linear  stiffness  matrix 
[Kt]  =  nonlinear  stiffness  matrix 
[Ka]  =  initial  stress  matrix 

The  full  expression  for  the  tangent  stiffness  matrix 
can  now  be  written  as: 

[Kt]  =  [Ko]  +  [Kl]  +  [Ka]  (2.31) 


Uith  an  expression  for  the  elements  contained  in  the  tangent 

®  stiffness  matrix  (Eqn  (2.31)),  the  derivation  can  proceed  in 

determining  the  terms  in  each  of  the  matrices  that  make  up 
the  tangent  stiffness  matrix. 

•  After  formulating  the  kinematic  relation  incorporating 
linear  and  nonlinear  in-plane  and  bending  strains  (see 
Appendix  F,  Eqn  (F.ll))  [37],  introducing  the  definition  for 

•  the  material  matrix  (see  Appendix  F,  Eqn  (F.12))  [37],  and 
determining  the  displacement  relations,  the  expression  for 
the  linear  and  nonlinear  stiffness  matrices  (  [Ko]  and  [Kl]  ) 

•  are  [37] : 


-  J, 
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and 


(zKri)  =  L 


[1]  [2] 
[2]  [3] 


dA 


(2.33) 


where , 

[b'I^Ca]  [Bj*]  +  CBf3^[A]  [bP]  +  [b|1^[a]  [bP] 

(12x12)  o  L  L  L  L  o 

[BP]^CA][Bf]  +  [Bf3^CA]CBfl  +  [B|1^[B][B5  (2.34) 

(l2xl2)  o  L  L  L  L  0 

[Bh  +  [Bh‘'[B]  [b;1  +  [b;^^[B]  [B**] 

(I2xl2)  o  L  L  L  L  o 


•  The  final  expression  necessary  for  the  tangent 
stiffness  matrix  is  the  initial  stress  matrix  [Ka]  .  The 
initial  stress  matrix  can  be  broken  down  into  in-plane  and 

•  bending  components.  Each  of  these  components  are  then 
defined  in  terms  of  the  appropriate  stress  resultants  and 
displacement  relations  (see  Appendix  F,  Eqns  (F.35)  to 

•  (F.42))  [37].  The  final  expressions  for  each  initial  stress 
matrix  component  are  (Appendix  F,  Eqns  (F.43)  and  (F.44)) 
[37]  : 


N 

xy 

N 

xy 

0  N 


0 

0 


N 


0 

0 

N 

I 

N 


xy 


dA 


(2.35) 


xy- 


2-25 


and 


(2.36) 


Uith  Eqns  (2.35)  and  (2.36),  the  full  expression  for 
the  initial  stress  matrix  can  be  written  as: 


( 


2.2.5  Plate  Elements  Representing  a  Shell 


(2.37) 


As  stated  earlier,  STAGS,  uses  plate  elements  to 
represent  the  surface  of  the  shell.  The  quadrilateral  plate 
elements  primarily  used  in  this  thesis  are  the  STAGS  QUAF 
411  elements.  The  QUAF  411  element  has  32  degrees  of 
freedom;  four  rotations  and  three  displacements  at  each 
corner  node  and  an  in-plane  displacement  at  each  of  the  four 
side  nodes.  This  element  is  shown  in  Fig  2.5. 

From  examining  Fig  2.5  it  is  noted  that  a  degree  of 
freedom  not  usually  used  in  plate  elements  is  included:  the 
normal  rotation  0  .  This  normal  rotation,  for  the  QUAF  411, 

Z 

is  the  average  in-plane  rotation  of  the  two  adjacent  edges 
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of  each  element  are  transformed  into  the  same  reference 
system,  compatibility  is  given  as  [42]: 


[0'  -  0^ 
1  y  yJ 

cos(  a/2  )  - 

[0’  +  0"1 
L  * 

[<  - 

cos(a/2)  - 

sin(  a/2  ) 

sin(  o/2  ) 


0 


0 


(2.38) 


where  the  superscripts  1  and  2  are  the  associated  element 
numbers.  As  the  angle  between  the  elements,  a,  becomes 
smaller  the  system  of  equations  in  Eqn  (2.38)  becomes 
increasingly  ill-conditioned  and  is  singular  at  a  -  0  [42] . 
Once  this  limit  is  reached,  6  is  omitted  and  Eqn  (2.38)  is 

Z 

assumed  as  0^  =  0^.  The  system  of  equations  becomes  flat 
y  y 

plate  relations. 

Another  problem  associated  with  the  use  of  flat 
elements  representing  a  curved  surface  has  to  do  with 
displacement  conformity  at  an  interface  of  two  flat 
elements.  For  a  flat  element  the  lateral  deflection,  w,  is 
usually  represented  by  at  least  a  cubic  polynomial  in  order 
to  handle  the  second  derivatives  associated  with  bending 
strain.  The  in-plane  displacements,  u  and  v,  are  usually 
represented  with  a  linear  or  quadratic  polynomial.  If  the 
compatibility  relation  for  two  adjacent  flat  elements  is 
derived,  the  resulting  expression  is  given  as  [42]: 
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Fig  2.6  Flat  Plate  Element  Representing  a  Curved  Shell  [42] 


v^j  cos(oe/2) 
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(2.39) 


If  one  examines  Eqn  (2.39)  it  is  evident  that  the 
displacement  compatibility  along  the  interface  of  the  two 
elements  cannot  be  satisfied  (if  v  is  quadratic,  it  cannot 
fit  the  same  curve  as  the  cubic,  m)  unless  v  and  w  are 
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represented  by  polynomials  of  the  same  order.  To  overcome 
this  displacement  nonconformity  the  QUAF  411  uses  a  third 
order  polynomial  to  represent  u  in  the  y-direction 
(quadratic  in  y)  and  v  in  the  x-direction  (quadratic  in  x). 
The  lateral  deflection,  w,  is  represented  by  a  cubic 
polynomial  in  x  and  y  .  The  QUAF  411  element  adds  a  shear 
term,  7,  at  each  corner  node  and  an  in-plane  tangent 
displacement,  t,  at  the  side  nodes  [42].  Rotational 
compatibility  is  enforced  only  at  the  nodes,  thus  the  411 
element  is  a  nonconforming  bending  element.  The  derivation 
of  the  shape  functions  associated  with  these  assumed 
displacement  fields  is  given  in  [42] . 

2.2.6  Equations  of  Motion 

Since  the  analysis  is  static,  the  problem  is  one  of 
solving  the  equation  of  static  equilibrium,  namely: 

^  Forces  *  0  (2.40) 

During  static  equilibrium,  the  potential  energy,  n  ,  must  be 

p 

at  a  minimum  and  the  first  variation  of  potential  energy, 

5n  ,  will  be  zero.  This  yields  the  equilibrium  equation. 

p 

For  bifurcation  analysis  (either  linear  or  nonlinear),  the 
second  variation,  5^n  ,  is  also  zero  at  the  collapse  point. 

p 
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The  nonlinear  studies  considered  in  this  effort  seek  a  limit 
point  in  the  loading  curve  rather  than  an  actual  bifurcation 
point.  It  is  well  known  that  curved  cylindrical  panels  have 
post  buckling  at  decreased  loads  only,  so  the  first  limit 
point  is  the  actual  collapse  load,  as  seen  in  Fig  2.7.  As 
shown  in  Fig  2.7,  the  solution  derived  by  linear  bifurcation 
may  be  either  higher  or  lower  than  the  actual  solution  [44] . 
For  a  cylindrical  panel  with  cutouts,  the  bifurcation  value 
is  too  high.  While  with  large  cutouts  the  bifurcation 
values  become  too  low  [35]  .  The  use  of  nonlinear 
bifurcation  is  limited  and  will  not  be  used  in  this  thesis. 

The  total  potential  energy  in  a  body  is  the  internal 
strain  energy  minus  the  work  of  the  applied  forces,  or: 

n  =  U  -  W  (2.41) 

p 

Both  U  and  U  can  be  calculated  based  upon  the 
displacement  field,  with  (Eqn  (C.5),  Appendix  C): 


U  = 


1 

2 


e®)^[D](e®)  dA 


(2.42) 


where  (e**)  and  [D]  are  defined  in  Appendix  D  and  reference 
[36].  The  work  done  by  the  applied  forces  is  (Eqn  (C.12), 
Appendix  C): 
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Displacement 


Fig  2.7  Load  Displacement  Curve  for  Typical 
Cylindrical  Panel 
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u  = 


(a)^(f  ) 


(2.43) 


where , 


I 


(a)  =  vector  of  displacements 
(f)  =  vector  of  externally  applied  forces 
For  applied  displacements.  STAGS  calculates  a  load  which 
results  in  the  proper  displacements. 

The  details  of  the  modified  Newton-Raphson  method  are 
given  in  the  STAGS  theory  manual  [26] ,  by  which  a  stiffness 
matrix  from  a  previous  step  is  used  to  estimate  the 
displacement  for  a  given  load  increment. 

Although  Bauld  [16]  used  a  finite  difference 
formulation  in  solving  the  strain  energy  relationship,  the 
method  in  determining  the  equation  of  motion  is  conceptually 
the  same  as  STAGS.  Bauld  separated  the  strain  energy  into 
quadratic,  cubic,  and  quartic  displacement  terms: 


2U  =  daa  -i-d  aaa  -t-d  aaaa  (2.44) 

Ij  i  j  ijk  i  j  k  Ijkl  1  J  k  1 

The  tensors  d  ,  d  ,  and  d  are  independent  of  the 
displacement  vector  a^  and  are  assumed  to  be  symmetrized. 
Details  for  the  tensor  determination  are  found  in  [16] . 

Bauld  goes  on  to  give  the  total  potential  energy  as: 
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n  =  ft  A 

P  ^2  ri 


+  i-  N1 

o  rf 


+  ^  N2  a  a 

IZ  r  8  i  r  i 


-  R  a 

s  a 


(2.45) 


A  is  the  linear  stiffness  matrix  with  no  dependence  on  the 

r  8 

displacement  vector,  a.  N1  and  N2  are  stiffness 

rs  rs 

matrices  with  linear  and  quadratic  dependence,  respectively, 

on  displacement.  R  is  the  applied  surface  force  vector. 

8 

Taking  the  first  variation  of  this  in  order  to  satisfy 
minimum  potential  energy,  one  obtains  the  equilibrium 
equation: 


Ta  +  4  N1  +  4  N2  la  -  R  =0  (2.46) 

ra  2  ra  J  raj  r  a 

2.2.7  Solution  Techniques 

The  STAGSC-1  computer  program  uses  a  Newton-Raphson  or 
a  modified  Newton-Raphson  (user  defined)  solution  technique 
to  solve  the  nonlinear  equilibrium  equations.  For  a  linear 
problem  in  STAGSC-1 ,  the  solution  technique  is  much  easier 
since  the  stiffness  matrix  is  not  a  function  of 
displacements.  The  linear  set  of  equations  is  given  as: 


[K](a)  =  (f) 


(2.47) 
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where  [K]  is  the  stiffness  matrix  (constant),  (a)  is  the  set 
of  nodal  degrees  of  freedom  (nodal  parameters),  and  (f)  is 
the  applied  loads.  STAGSC-1  uses  a  Cholesky  triangular 
decomposition  with  forward  and  backward  substitution  to 
solve  Eqn  (2.47)  for  (a)  [42].  The  solution  for  the 
nonlinear  equilibrium  equations  is  much  more  involved  and 
is  considered  next. 

An  explanation  of  the  reference  systems  in  STAGSC-1  is 

necessary  before  progressing.  There  are  basically  two 
reference  systems  defined:  a  global  reference  system 
(subscript  g)  and  a  local  element  reference  system 
(subscript  e).  The  global  reference  system  is  fixed  in 
space  and  does  not  move.  The  local  reference  system  uses 
and  updated  Lagrangian  approach.  The  local  reference 
system,  called  a  corotational  system,  is  fixed  to  the 
element  and  moves  with  the  element  during  rigid  body  motion 
[25] .  This  allows  for  the  removal  of  rigid  body  motion  in 
the  element  before  calculating  strain.  The  present  version 
of  STAGSC-1  (1986)  also  redefines  the  standard  way  of 
representing  a  rotation  as  a  vector  quantity.  Rotations  are 
defined  by  a  triad  (three  mutually  perpendicular  unit 
vectors)  to  accurately  map  local  rotations  since  they  depend 
on  order  of  rotation.  The  previous  versions  of  STAGS  used 
vectors  to  describe  rotations  which  with  large  rotations 
gave  spurious  results  (vectors  cannot  track  order  of 


rotations  and  are  only  good  for  small  rotations).  This  is 
the  reason  the  QUAF  411  element  was  developed.  The  411 
element  handles  larger  rotations,  but  with  increased  cost 
[46] .  The  reader  is  referred  to  [47]  for  a  complete 
discussion  on  this  rotational  formulation.  Finally, 
differentiation  and  integration  within  the  element  are  done 
with  respect  to  this  corotational  system  [25] . 

For  the  STAGSC-1  quadrilateral  element,  the  local 
reference  system  is  defined  as  follows:  An  approximation  of 
the  warped  element  surface  is  made  by  crossing  the  principal 
diagonals  of  the  element  to  form  a  vector  and  the 
establishing  a  plane  normal  to  this  vector  such  that  one  of 
the  nodes  lies  on  this  plane.  This  node  is  where  the  local 
reference  system  is  established  by  taking  the  z-axis  normal 
to  the  plane,  the  x  or  y  axis  is  located  along  one  of  the 
element  edges,  and  the  remaining  axis  completes  a  Cartesian 
right-handed  system  (see  Fig  2.8)  [47]  .  All  local 
deformations  of  the  element  are  referenced  to  this  local 
system  during  an  iteration  for  a  solution.  Once  a  solution 
has  been  found  in  the  global  system  the  local  coordinates 
are  updated  and  a  new  local  reference  system  is  established. 
This  movement  of  the  local  reference  system  appears  to  be 
Euler ian  except  that  the  local  coordinates  of  a  point  change 
[25] .  Strains  and  rotations  in  e  local  reference  system 
are  usually  small  enough  so  as  to  ignore  the  nonlinear 
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stiffness  matrix  [Kl]  and  sometimes  even  the  initial  stress 
matrix  [Ka]  that  are  a  part  of  the  tangential  stiffness 
matrix  [25].  This  is  only  true  for  the  local  system,  the 
global  tangent  stiffness  matrix  must,  in  general,  contain 
all  the  terms  (see  Eqn  (2.31)).  From  thorough  research  of 
the  STAGS  documentation,  it  appears  that  STAGS  uses  the  full 
tangent  stiffness  matrix  in  both  the  global  and  the  local 
systems  to  overcome  difficulties  in  solutions  due  to  highly 
nonlinear  problems. 
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With  the  reference  systems  defined,  the  Newton-Raphson 
solution  techniques  can  be  explained.  To  start,  expand  the 
expression  for  the  sum  of  the  external  and  internal  forces, 
(tp)  (see  Appendix  C),  in  a  Taylor  expansion: 


=  IvCCa)")! 


d(  v>) 
d(a)J 


(Aa)"  +  H.O.T.  =  0  (2.48) 


where  H.O.T.  means  higher  order  terms  and, 

(a)"*^  =  (a)"  +  (Aa)"  (2.49) 

From  Eqn  (F.4)  (see  Appendix  F): 

=  [Kt]  (2.50) 

d(  a  ) 

where  [Kx]  is  the  tangent  stiffness  matrix,  a  function  of 
displacements.  If  one  uses  the  expression  for  (yj)  in  Eqn 
(C.15)  (Appendix  C)  and  writes  it  as  a  function  of 
displacements  for  the  nonlinear  problem  one  obtains  [27] : 

(V<a))  =  [K(a)](a)  -  (f)  (2.51) 
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where 


[K(a)](a)  =  internal  resisting  forces  of  the  element 
(f)  =  externally  applied  loads 
Inserting  Eqns  (2.50)  and  (2.51)  into  Eqn  (2.482)  and 
eliminating  the  higher  order  terms  results  in: 

[KT](Aa)''  =  (f)  -  [K(a)"](a)"  (2.52) 

The  problem  now  is  to  find  the  displacements  within  the 
elements  that  balance  the  externally  applied  loads  and  the 
internal  resisting  forces  (i.e.  the  sum  of  the  right  hand 
side  of  Eqn  (2.52)  equaling  zero).  Since  an  updated 
Lagrangian  approach  is  being  used,  the  internal  deformations 
are  relative  to  the  local  reference  system,  therefore  the 
last  expression  in  Eqn  (2.52)  needs  to  be  written  in  terms 
of  the  local  displacements.  Rewriting  Eqn  (2.52)  in  terms 
of  the  local  (subscript  e)  and  global  (subscript  g) 
displacements  results  in  [25] : 

[KT](Aa'’)  =  (f)  -  y  [k(a  )"](a  )"  (2.53) 

9  4,  e  o 


where , 


^  [k(a^)''](aj" 


=  internal  resisting  forces  transformed 
into  the  global  system  and  summed 


For  simplicity,  the  Newton-Raphson  solution  technique  is 
given  as  follows  [25] : 

1.  Increment  the  applied  load  (f). 

2.  Establish  the  local  coordinates  for  each  element. 

3.  Formulate  the  element  tangent  stiffness  matrices 
in  terms  of  their  local  degrees  of  freedom. 

4.  Transform  the  local  tangent  stiffness  matrices  to 
the  global  coordinate  system  and  assemble  them  into  the 
global  stiffness  matrix. 

5.  Compute  the  values  of  the  local  degrees  of  freedom 
(a  )  (zero  for  the  first  iteration  within  each  load 

e 

step )  from  the  global  degrees  of  freedom  (a  ) . 

g 

6.  Calculate  the  element  internal  forces  [k(a  )](a  ). 

e  a 

Transform  these  forces  to  the  global  system  and 

assemble  them  with  the  other  element  internal  forces. 

7-  Solve  Eqn  (2.53)  for  ( Aa  ). 

g 

8.  If  the  vector  ( Aa  )  is  not  small  enough  (i.e. 

g 

converges)  return  to  step  3. 

After  the  solution  converges  (step  8)  return  to  step  1  and 
repeat.  Once  the  solution  has  converged  for  the  final  load 
step  ((f)  has  been  incremented  to  equal  the  final  load)  the 
solution  is  complete. 

An  alternate  solution  technique  is  the  modified 
Newton-Raphson  method.  This  is  the  method  used  in  STAGS  for 
this  effort.  In  this  method  the  tangent  stiffness  matrix  is 


not  reformulated  at  every  iteration  within  a  load  step.  The 
previously  assembled  stiffness  matrix  is  used  for  successive 
iterations  and  is  only  refactored  when  the  convergence  rate 
dictates  [25] . 

STAGS  allows  the  user  to  control  many  of  the  parameters 
in  the  solution  process  to  include:  either  using  the  full 
or  modified  Newton-Raphson  method,  control  of  the  load  step 
size,  the  number  of  attempts  at  a  solution  if  convergence 
difficulties  arise,  and  whether  the  Riks  solution  algorithm 
or  the  Thurston  Equivalence  Transformation  method  is  used. 
For  this  effort,  the  Riks  solution  algorithm  is  used.  This 
algorithm  is  explained  in  greater  detail  in  Section  2.3.5. 
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2.3  SHELL  Theory 


SHELL  was  developed  in  1986  by  Scott  Dennis  as  part  of 
his  Ph.D.  dissertation  [41] .  Its  primary  purpose  is  the 
analysis  of  thin,  cylindrical,  shell  structures.  It  is 
based  entirely  upon  the  finite  element  method.  SHELL’S 
primary  capability  is  the  nonlinear  large  displacement  and 
moderate  rotation  characteristics  of  thin  shells  with  the 
added  feature  of  including  through  the  thickness  shear 
strain.  The  basic  element  used  in  SHELL  is  a  cylindrical 
shell  element  that  can  exactly  match  the  curvature  of  the 
shell  surface.  The  nonlinear  analysis  is  based  upon  either 
the  Donnell  approximations  (similar  to  Von  Karman  plate 
relations  but  has  curvature  incorporated  within  the 
kinematic  equations)  or  the  exact  nonlinear  kinematic 
relations.  And  the  additional  feature,  mentioned  above,  of 
this  program  is  that  SHELL  incorporates  a  parabolic 
distribution  for  transverse  shear  strain. 

For  this  effort,  the  exact  nonlinear  kinematic 
relations  are  discussed  herein  because  of  the  larger 
rotations  (on  the  order  of  e  30®)  encountered  in  the  panels 
that  are  not  simply  supported  along  the  vertical  edges  (see 
Sections  3  and  4).  These  large  rotations  far  exceed  the 
assumptions  made  for  the  Donnell  relations,  thereby 
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invalidating  them.  A  brief  look  at  the  theory  for  this 
cylindrical  shell  element  and  the  governing  relations  used 
in  SHELL  is  presented  in  the  following  sections. 

Bathe  and  Ho  [48]  stated  the  following  set  of  criteria 
for  a  desirable  shell  element: 

1 .  )  No  spurious  zero-energy  modes  should  exist,  so  that 

reliable  results  can  always  be  expected.  No 
numerical  fudge  factors  should  be  necessary, 
either . 

2. )  The  element  should  be  applicable  to  general  shell 

structures,  including  those  with  beam  stiffeners, 
cutouts,  intersections,  etc. 

3. )  The  element  should  be  cost-effective  for  linear  as 

well  as  nonlinear  static  and  dynamic  analysis. 

This  implies  that  the  degrees  of  freedom  is  held  to 
a  minimum.  It  should  allow  analysis  of  large 
displacement  and  large  rotation  problems,  and 
materially  nonlinear  situations. 

Criterion  1  is  achieved  in  this  formulation.  The  code  has 
not  been  developed  beyond  the  research  stage,  so  application 
to  other  than  plate  and  cylindrical  shell  cases  modeled  with 
rectangular  elements  is  not  yet  possible .  However ,  the 
theory  is  not  restricted  geometrically,  so  criterion  2  is 
somewhat  met.  Material  linearity  is  assumed  in  the 
formulation  of  the  element,  so  criterion  3  is  only  partially 
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fulfilled.  The  incorporation  of  through  thickness  shear 
effects  in  a  shell  structure  while  maintaining  a  two 
dimensional  analysis  well  satisfies  the  first  part  of  this 
criterion,  however. 

The  following  development  applies  specifically  to 
cylindrical  geometry.  Dennis  [24]  presents  the  general  case 
of  a  doubly  curved  shell. 

2.3.1  Geometry  and  Assumptions 

The  curvilinear  orthogonal  coordinate  system  and 
nomenclature  used  in  this  formulation  of  the  laminated 
cylindrical  shell  is  shown  in  Fig  2.9.  The  x-axis  lies 
along  the  straight  dimension  of  the  panel;  the  s-axis 
follows  the  circumference,  and  the  z-axis  is  everywhere 
normal  to  the  shell  middle  surface,  positive  toward  the 
center  of  curvature.  The  surface  formed  by  the  x  and  s  axes 
lies  in  the  center  of  the  thickness  of  the  panel,  so  the 
thickness  coordinate  is  negative  on  the  outer  surface  and 
positive  on  the  inner  surface.  Displacements  along  the  x,  s 
and  z  axes  are  u,  v  and  w  respectively.  In  the  early 
vectorial  development  the  coordinate  Q  is  used  instead  of  s 
for  the  purpose  of  generality,  s  is  used  when  specializing 
to  the  cylindrical  geometry.  Since  the  structure  analyzed 
here  is  an  open  shell,  the  angle  8  is  also  useful  for 


describing  shallowness.  The  angle  <|>  specifies  the 
orientation  angle  of  each  ply  in  the  laminate.  Subscripts 
denoting  stress  and  strain  orientation  are  explained  in 
Table  2.2  and  Figure  2.4. 

Table  2.2.  SHELL  Contracted  Notation 


stress 

strain 

cylindrical 

explicit Icontracted 

explicit Icontracted 

coordinates 

a  a  ^ 

11  1 

e  e 

11  1 

X  -»  1 

a  a 

22  2 

e  e 

22  2 

s  2 

a  a 

33  3 

e  e 

33  3 

z  -»  3 

shear : 

a  a 

23  4 

e  e 

23  4 

s-z  -♦  4 

a  a 

13  5 

13  S 

x-z  -♦  5 

a  o 

12  6 

e  e 

12  6 

x-s  -4  6 

2.3.2  Constitutive  Development 

The  development  of  the  constitutive  relations  for  SHELL 
is  shown  in  Appendix  E.  It  is  identical  to  Section  2.2.2 
except  it  assumes  a  modified  plane  stress  relationship. 

This  modified  relationship  allows  a  =  0  but  it  assumes 
and  are  not  equal  to  zero,  thereby  incorporating  nonzero 
through  the  thickness  shear  stress.  Therefore,  following 
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the  same  logical  progression  in  Section  2.2.2,  but  now 
incorporating  the  modified  plane  stress  relation,  the 
reduced  stiffness  constitutive  relation  becomes  [45] : 


Q  Q  0  0 

11  12 

Q  Q  0  0 

12  22 

0  0  Q  0 

66 


(2.54) 


Q  0 

4  4 


where  the  Q  are  reduced  stiffness  coefficients  related  to 
i  J 

the  S’s  by 


S  S 

Q  =  S  - 


(2.55) 


In  terms  of  the  engineering  coefficients  again. 


V  E 
21  2 


(2.56) 
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Mhere,  A  =  1  -  u  u  . 

12  21 

Finally,  in  order  to  analyze  a  stack  of  plies,  they 
must  all  be  referenced  to  a  global  axis  system  and  their 
effects  summed  (Appendix  E,  Eqn  (E.6)): 

R=  [  ^  ]  [°u],  [  ^  TR 

Mhere  (see  Appendix  E,  Eqn  (E.5)), 


and  c  =  cos(d),  8  =  sin(9).  Uith  this  transformation  the 
constitutive  relations  are: 


and  the  transformed  reduced  stiffnesses  are  shown  in 
Appendix  E,  Eqn  (E.8). 


2.3.3  Strain-Displacement  Relations 

The  geometric  nonlinearity  arising  from  the  panel’s 
curvature  is  incorporated  by  the  strain-displacement 
relations.  It  is  also  here  that  the  through  thickness  shear 
effects  are  incorporated  into  the  analysis.  This 
development  follows  Dennis  [41]  who  authored  the  computer 
code  which  employs  these  relations. 

The  incorporation  of  transverse  shear  effects  is 
accomplished  by  assuming  a  modified  state  of  plane  stress 
for  the  lamina  in  which  a  =  O  (and  hence  c  =0)  but  a  and 

a  3  4 

a  are  allowed  to  be  small  nonzero  values.  These  transverse 
s 

shear  stresses  are  assumed  to  equal  zero  on  the  top  and 
bottom  surfaces  of  each  ply.  and  the  associated  strains  will 
vary  parabolically  through  the  ply  thickness. 

To  proceed,  keeping  Just  the  linear  (first  order) 
displacement  terms  for  the  transverse  shear  strains  and 
c  we  see  that  [2] : 

5 

*■4  ~  R”  [^3*2  *  ^2^2*3  ^2^2*3] 

^  (2.60) 

e  "c-u*  +hu,  -uh 

5  FT  3’i  1  1*3  11 

where  the  h^  are  the  coordinate  system  scale  factors;  for 
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1  and 


the  cylindrical  geometry  used  in  this  work,  h^  = 

h  =  1  -  z/R. 

2 

The  displacement  equations  in  the  thickness  variable  z 
which  permit  the  incorporation  of  the  desired  through 
thickness  feature  are  [41] : 

u(x,d,z)  =  u°  +  zv^  + 

v(x,0,z)  =  v®^l  -  +  zt|>2  +  z^a^  +  +  z^e^  (2.61) 

w(  X  ,  0  )  =  w 

where  u",  v“,  w,  ,  0^,  and  0^  are  functions  of  the 
coordinates  x  and  0.  The  displacements  u®  and  v°  are  of  the 
shell  middle  surface;  transverse  displacement  w  is  the  same 
throughout  the  thickness  since  transverse  normal  strain  is 
assumed  negligible.  The  ip^terms  are  rotations  of  the 
surface  normals  in  the  x  and  8  planes,  and  0  ,  y  and  9  are 

i  i  1 

to  be  found  by  applying  the  assumption  that  transverse  shear 
stresses  and  are  zero  on  the  shell  surfaces. 
Substituting  the  equations  for  v  and  w  into  Eq  (2.60)  for 
yields  [41] : 


^4  1-z/R 


f  ® 

+  ( 1-Z/R)  1^-  ^  +  v>2  +  2z02  +  az^y^  +  4z^0. 


1  -  I  +  ZV^  -H  -H  Z^^  -K  z\ 


(2.62) 


For  zero  transverse  shear  stress  at  the  surfaces,  the 
associated  strain  will  also  be  zero.  Enforcing  this 
condition  by  substituting  ih/2  for  z  in  Eqn  (2.62),  setting 
both  resulting  expressions  to  zero  and  hence  equal  to  each 
other,  then  solving  for  the  unknown  variables  produces  [41] 


0 


=  0 


e.  =  ^ 


2 

2R 


(2.63) 


1  - 


8R‘ 


-4 

3h' 


(»2  *  “’2] 


Note  that  an  h/R  value  of  1/5  (quite  high  for  practical 

aerospace  shell  geometries)  allows  neglect  of  the  h/R  term 

in  the  left  side  of  Eqn  (2.63).  Replacing  0  ,  0  and  y  in 

22  2 

Eqn  (2.62)  with  the  results  in  Eqn  (2.63)  yields  the 

transverse  shear  strain  e  in  terms  of  transverse 

4 

displacement  w  and  rotation  [41]  : 


^“-2  ^  »2> 


1  - 


4z 


8z 


3h^R 


(2.64) 
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Note  again  that  there  are  1/R  terms  to  be  neglected. 
Therefore,  the  transverse  shear  strain-displacement 
relations,  e  ,  along  with  a  similar  analysis  with  e 

4  5 

(simpler  since  h^  =  1),  become  [41]: 

(2.65a) 

(2.65b) 

Replacing  0^,  0^  and  in  Eq  (2.61)  yields  the 
displacement  equations  [41]  : 


^4  1-z/R 


<w,2  + 


1  - 


4z‘ 


=  (w,  +  tt»  ) 

D  11 


1  - 


4z 


u(x,6,z)  =  u®  +  zm  -  z^ 

‘  3h^ 

v(x.0,z)  -  v’[l  -  ^]  t  z»^  -  ^  z=[»^  (2.66) 

w(  X ,  e )  =  w 

At  this  point  one  can  note  that  this  formulation 
provides  seven  degrees  of  freedom:  u'  v;  w;  w,^;  w,^; 
and  w  . 

Now  that  we  have  displacement  equations  which 
incorporate  i  parabolic  through  thickness  shear  stress 


*  "-i] 
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distribution,  the  in-plane  kinematic  equations  are  derived 
for  the  shell  middle  surface.  The  fully  general  strain 
displacement  relations  are  quite  extensive;  the  theoretical 
development  incorporating  them  can  be  found  in  Dennis’  large 
displacement  -  moderate  rotation  general  shell  development 
[41] .  The  general  strain  displacement  relations  of  Eqns 
(2.9)  and  (2.10)  with  the  kinematics  of  Eqn  (2.66)  will  give 
the  in-plane  shell  strain  displacement  relations.  These 
expressions  can  be  specialized  for  a  shell  geometry  of 
interest  by  defining  the  scale  factors,  h  ,  used  in  Eqn 

i 

(2.67)  [41]  . 

(2.67) 

2 

where  a  =  g  (no  sum ) . 

y  yy 

Substituting  Eqn  (2.67)  into  Eqns  (2.9)  and  (2.10)  and 


h  =  a  ( 1  -  z/R  ) 
11  1 

^  -  0/1  -  ^/R3) 


carrying  out  the  indicated  differentiation  yields  for  the 
in-plane  strain-displacement  equations  [41] : 


e°  +  zx^  +  z^x^  +  z^x^  +  z^x^  +  z®x® 


c°  +  zx'  +  z^x^  +  z\^  +  z^x^  +  z®x® 

2  2  2  2  2  2 


(2.68) 


h  h 
1  2 


1^  22^  33_^  4  4  66 

e  +ZX  +  ZX,  +  ZX  +  ZX  +  ZX 
6  6  6  6  6  6 


where  the  e°  and  the  xj  terms  (j«i,2,6;  i«i,2,3,4,6)  are 
functions  of  the  displacements  and  the  scale  factors,  and 
can  be  found  in  Appendix  A  of  [41] .  It  is  noted  that  the 
superscripts  on  the  xj  are  not  exponents,  they  are 
individual  strain  components  that  correspond  to  the  power  of 
z  that  multiplies  it  ,  and  the  subscripts  on  x^  indicate  the 

strain,  e  ,  e  ,  or  e  that  these  components  correspond  with. 

12  6 

The  following  equivalent  representation  of  the  strains 
is  conducive  to  the  matrix  operations  of  the  potential 
energy  formulation  in  the  next  section.  The  in-plane  strain 


displacement  relations  are  represented  by  [41] : 
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(2.69a) 


where , 


'1  “  ^  *^’0 


where , 


S  =  '"2 


e  I  +  ;  (p  .  1....7) 


g  +  ifu?  +  ^  +  !i!1 

R  2[  2  2  2  ^2  ^2) 


(2.69b) 


VW  ,  V  ,  w 

2  2 


e  =  e®  +  z‘’x  ;  (p  -  i...,7) 

6  6  6p 


(2.69c) 


where , 

e°  =  u,  +  v,  +  u,u,  + 
6  2  1  12 


V,  V,  w,  w,  +  ifvw,  -  V,  w] 

’l  ’2  ’l  ’2  Rl  1  1  J 


Note:  The  higher  order  bending  strain-displacement 

relations,  x  (1-1,2, 6;  p-i,..,?),  are  completely  shown  in 

ip 

[41]  . 


The  transverse  shear  strains  (in  similar  form)  are 


[41]  : 


O  D 

e  =  e  +  z*^x^  ; 

44  4p 


e  =  w ,  +  IP 

4  2  ^2 


=  3k(w.^  + 


(2.70a) 


X  ( p  =  1 , 3, 4, 5, 6, 7) 

<p 


e  =  e®  +  z‘’x  ; 

55  5p 


e  =  w ,  +  u; 

5  1  ^1 


Xg2  =  3k(w,^  +  v^)  (2.70b) 


x_  (p  =  1,3, 4, 5, 6, 7) 

Sp 
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This  allows  assembly  of  the  strain  displacement  equations 
into  a  matrix  format  which  will  become  convenient  in  the 
next  section: 


£ 

1 

1 

'x 

1 1 

X 

12 

X 

13 

X 

14 

X 

15 

X 

16 

X 

17 

£ 

2 

►  S  1 

o 

S 

.  + 

X 

21 

X 

22 

X 

23 

^*24 

X 

25 

X 

26 

X 

27 

£ 

1  eJ 

o 

X 

^  61 

X 

62 

X 

63 

^64 

X 

65 

X 

66 

X 

67-* 

(2.71) 


In  a  general  expression, 

in  =  {8“)  +  [X  ]{£  ) 


(2.72) 


It  is  worth  noting  that  these  kinematics  avoid  the  common 
pitfall  of  shear  locking,  wherein  the  model  becomes  artificially 
stiff  as  the  shell  thickness  is  decreased.  This  is  a  problem 
with  finite  element  formulations  which  incorporate  constant  or 
linearly  distributed  shear  strain  through  the  thickness,  and  it 
necessitates  the  use  of  a  correction  factor .  However , 
examination  of  the  compatibility  relations  associated  with  the 
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strain  displacement  equations  developed  in  this  section  shows 
that  the  terms  associated  with  transverse  shear  drop  out  as 
thickness  is  reduced  to  zero  (see  [41] ). 

Note  the  kinematics  developed  in  this  section  are 
specifically  for  the  cylindrical  geometry.  The  analysis 
accomplished  for  this  thesis  uses  the  large  displacement  - 
moderate  rotation  formulations. 

2.3.4  Equations  of  Motion 

The  shell  potential  energy  is  the  sum  of  the  internal 
strain  energy  and  the  work  done  by  external  forces: 

n  =  U  +  V  (2.73) 

p 

where  the  internal  strain  energy  is  given  by: 

U  =  do  dh  (2.74) 

where  Q  represents  the  shell  middle  surface.  The  internal 
strain  energy  U  is  composed  of  in  plane  and  normal  terms 
(set  1),  and  transverse  shear  terms  (set  2).  Expanding  the 
expression  for  and  by  inserting  Eq  (2.69)  through 
(2.72)  into  Eq  (2.74)  gives 


u 

1 


u 

2 


^  23, ^Ce”  +  *  3^^(£°  +  z-x^^y 


+  2Q  (  e®  +  z‘’x  )(  e®  +  z/x  )1 

262  2p  6 


dz  dQ 


(2.75) 


^  r  f  [q  (e®  +  z^x  )'"  +  Q  (  e°  +  z‘‘x  )‘ 
2  J_J  1  44  4  42'  55  5  52 


'Q‘'h 

t  20,^(£;  »  t  z'x^^)] 


dz  do 


where  p,r  =  1,2,..., 7.  Integrating  the  z  over  ±  h/2  yields 
the  equation  for  strain  energy  as  a  function  of  the  middle 
surface  only,  which  is  the  desired  formulation  for  this 
shell  problem.  A  further  simplification  performed  in  the 
coding  of  the  SHELL  program  is  one  of  symmetry  in  ply  lay-up. 
This  results  in  the  cancellation  of  elasticity  arrays  which 
are  multiplied  by  odd  powers  of  the  transverse  coordinate  z. 
Rearranging  (see  [41])  yields  the  final  forms  of: 


U 

1 


I  r  (S®)^[A]{S®)  dQ 

Jq 

+  I  J^2{S®>^^[B]  +  [D]  +  [E]i-[F]  +  [G]  +  [H]  +  [I]j{X)  dQ 
+  I  J^{X)^^[D]  +  [E]  +  [F]  +  [G]  +  [H]  +  [I]  +  [J]-^[K]i-[L] 

+  [P]-*-[R]  +  [S]  +  [T]  j{X}  dQ  (1,2,6) 


(2  76a) 
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(2.76b) 


U2  =  I  J  [{8®>^[A]{X>  +  2{S°)[D][X] 
+  CX]’'[F][X]j  do 


(4.5) 


where  {[A,  B.  D,  E,  F,  G,  H,  I,  J,  K,  L,  P,  R,  S,  T]} 


1  2 


1 3 


z'Mdz. 


(2.77) 


2.3.5  Finite  Element  Formulation 

The  solution  of  the  static  finite  element  problem 
involves  finding  the  equilibrium  state  between  applied  load 
and  structural  response.  This  state  can  be  determined  by 
finding  where  the  system  potential  energy  is  stationary, 
i.e.  where  its  variation  is  zero.  Recall  Eq  (2.73): 

n  =  U  +  V  (2.78) 

p 

where  now  internal  strain  energy  can  be  represented  by 

U  =  I  a”^  [  *^2]  *  (2.79) 

where  a  -  column  array  of  nodal  displacements.  K  »  constant 
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stiffness  terms,  =  stiffness  terms  linear  in 
displacement,  and  =  stiffness  terms  nonlinear  in 
displacement.  This  collective  group  of  stiffness  terms  is 
known  as  the  tangent  stiffness  matrix  K  ,  since  it  defines 

T 

the  local  slope  of  the  load  -  displacement  curve  at  any 
point  of  equilibrium.  The  external  work  can  be  represented 
as 


V  =  -  a’^XP  (2.80) 

where  P  is  a  column  array  of  applied  nodal  loads  and  X  is  a 
multiplier  whose  utility  will  be  seen  in  the  section  on  the 
solution  algorithm.  Equilibrium  is  defined  as  the  state  at 
which  internal  strain  energy  and  external  work  balance; 
potential  energy  is  at  a  relative  minimum  here.  This  point 
can  be  found  by  substituting  Eqns  (2.79)  and  (2.80)  into  Eqn 
(2.82)  and  taking  the  first  variation: 

=  6a^  I^K^a  -  xpj  =  O  (2.81) 

Since  displacements  5a  are  nonzero  for  all  but  the  trivial 
solution,  the  bracketed  expression  must  be  zero  for 
equilibrium.  It  is  a  function  (call  it  g)  of  a  and  X: 

K  a  -  XP  =  0  =  g(a,X)  (2.82) 

T 
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since  K  varies  with  load  and  displacement,  a  numerical 

T 

iterative  solution  is  preferred  to  attempting  to  solve  in 
closed  form.  Therefore,  to  determine  the  nodal 
displacements  a  for  a  given  applied  load  P,  Eqn  (2.82)  is 
solved  iteratively  by  the  modified  Newton-Raphson  method  as 
discussed  next. 

The  technique  advanced  by  Riks  [49]  and  Mempner  [50] 
and  demonstrated  by  Crisfield  [55]  of  incrementing  a  desired 
arc  length  along  the  load-displacement  curve,  while  solving 
iteratively  via  the  modified  Newton-Raphson  method,  is  the 
solution  technique  used  in  this  work.  The  modified 
Riks-Wempner  algorithm,  added  to  the  SHELL  program  by  Tsai 
and  Palazotto  [51]  ,  allows  for  tracing  of  the 
load-displacement  response  through  both  load  reversing  (snap 
through)  and  displacement  reversing  (snap  back)  critical 
points,  so  the  most  complex  behavior  can  be  continuously 
followed. 

The  essence  of  the  Riks/Uempner  method  is  that  neither 
load  P  nor  displacement  a  is  independently  controlled; 
rather ,  a  selected  arc  length  of  the  load-displacement 
curve  is  incremented.  The  equilibrium  condition  is  found 
which  satisfies  the  constraint  relation: 


^  Aa  +  P^P  =  As^ 

i+l  l4t  1*1 


(2.83) 
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where  Aa  is  the  incremental  displacement  for  step  i+1 , 

i  4'! 

and  AX  is  the  fraction  of  load  P  applied  at  step  i+1  (see 

i4l 

Fig  2.10).  The  effect  of  constraint  Eqn  (2.83)  is  that  each 
subsequent  step  solution  is  searched  for  on  an  arc  of  radius 
As  from  the  current  solt.cion.  The  initial  value  for  the 
quantity  AX  is  specified  in  the  problem  data  deck. 

To  apply  the  modified  Newton-Raphson  method,  the  first 
variation  of  Eq  (2.84)  is  taken  and  applied  at  steps  i: 


K  5a  =  6X  P  -  g(a  ,X  ) 

T  i  i  i  i 


(2.84) 


where 


5a  =  5a  +  5X  5a  ^  (2.85a) 

i  i 1  1  i2 

and 

5a  =  -K"^g(a  ,X  )  ( out-of-balance  (2.85b) 

Til  term) 

5a  =  K'V  (linear  term)  (2.85c) 

i  2  T 


The  displacement  increment  snd  load  parameter 

increment  AX  are  updated  by 

i 


Aa  =  Aa  +5a 

i4l  i  i 

and 


AX  =  AX  +  5X, 

i  4l  1  i 

I 


(2.86a) 


(2.86b) 
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The  linear  solution  at  the  beginning  of  each  load  step 


is  defined  by 


Aa 

1 


AX  8a 

1  i2 


(2.87) 


where  6X^  is  the  incremental  loading  parameter  in  the 
beginning  of  each  load  step.  Crisfield  [55]  has  shown  that 
it  is  preferable  to  fix  the  “incremental  length"  As.  The 
constraint  relation  is  then  replaced  by 

Aa’'  Aa  ,  =  As^  (2.88) 

i4t 


for  all  iterations,  i.  By  substituting  Eqn  (2.86a)  into  Eqn 
(2.88),  one  obtains  the  following  quadratic  equation  for  the 


change ,  5X 


a5X  +  2b6X  +  c  =  0 

i  i 


• 

where 

a  =  6a 5a 

12  12 

b  =  ( Aa  +  6a  )’^5a 

1  11  1 

• 

c  =  (Aa  +6a^^)(Aa 

12 

1 


1 1 


-  As' 


(2.89) 
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The  total  displacement  a  and  the  total  load  parameter 

n 

\  for  the  nth  load  step  is  updated  by 

n 

a  =  a  +  Aa  ( 2 . 90a ) 

n  n- 1  n 

and 

X  =  X  +  X  (2.90b) 

n  1  n 

In  brief,  the  algorithm  proceeds  for  a  given  load  step  n  as 
follows  (i  denotes  the  iteration).  Refer  to  Fig  2.10= 

a. )  The  tangent  stiffness  matrix  at  the  current 

deformed  geometry  is  determined  and  held  constant 
throughout  the  iteration.  Select  a  reference  load 
P  and  reference  number  of  iteration  for  each  load 
step  as  N  ,  and  the  number  of  load  steps  M  .  A 

9  8 

load  parameter  X^  for  the  first  load  step,  a 
maximum  load  parameter  X  ,  and  a  convergence 

max 

criteria  c  for  each  load  step  are  also  selected. 

b. )  If  it  is  the  first  load  step: 

1)  Calculate  6a from  Eqn  (2.85c) 

2)  Calculate  the  constant  arc  length 

As  =  AX  (aa’’  6a  )’'^ 

1  0  12  12 

3)  Calculate  Aa^  from  Eqn  (2.87) 

4 )  Go  to  step  ( c-2 ) 


c.)  For  any  load  step  n: 

1)  For  the  first  iteration: 

(i)  As  =  As  (N  /  N  ) 

n  n-1  B  8-1 

(ii)  Calculate  5a^^  from  Eqn  (2.85c) 

(iii)  AX  =  ±  - - 

12  i  2 

AX^  is  positive  if  detCKx)  is  positive 
AX^  is  negative  if  det(KT)  is  negative 

(iv)  Calculate  Aa^  from  Eqn  (2.87) 

2)  For  any  iteration  i+1: 

(i)  Calculate  the  out-of-balance  term  6a 
from  Eqn  (2.85b) 

(ii)  Calculate  the  load  parameter  change  6X 


i  1 


from  scalar  quadratic  equation  (2.85). 

Two  roots  are  obtained  as  6X  and  6X 

i  1  i2 

If  6X.^  and  6X^2  are  imaginary,  let 
As  -As  /2,  and  go  to  step  (c-l-ii).  If 

n>l  n 

they  are  real  roots,  go  to  the  next  step. 

(iii)  Calculate  the  updated  displacement 
increment  Aa^  and  Aa^  from  Eqn 

14^1  i4l 

(2.86a)  for  the  two  roots  6X  and  6X 

il  12 

(iv)  Calculate  the  angle  0^  between  the  Aa^ 

and  Aa'  ,  and  9  between  the  Aa  and 
i4i  2  1 


Aa' 


i4l 

and  6X  =  6X 


If  9  is  positive,  Aa  =  Aa 

1  i+1  i+1 


i  1 


If  0^  is  positive. 
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i^a  .  Aa  and  6X  =  6X  .  If  both  9 

i  +  l  I'fl  i  i2  1 

and  0  are  positive,  6X  takes  the 

2  i 

value  closest  to  the  linear  solution 
5X=-c/2b . 

(v)  Update  the  load  parameter  AX  from  Eqn 

i 

(2.86b)  and  the  total  displacement 
a  =  a  +  Aa  -  Aa 

n  n  i  -f  1  i 

3)  Repeat  step  (2)  until  the  difference  betMeen 
the  norms  of  the  total  displacement  for  two 
consecutive  iterations  is  within  e  of  the 
norms  of  the  total  displacement.  The  number  of 
iterations  to  reach  the  convergent  criteria 

is  defined  as  N  ,  and  the  final  incremental 

n 

load  parameter  for  the  nth  load  step  is  AX  . 

in 

4)  Update  the  total  load  parameter  as  from  Eqn 
(2.90b). 

d . )  If  the  load  step  number  n  is  equal  to  M  or  the 

9 

total  load  parameter  X  is  larger  than  X  ,  stop. 

n  nax 

otherwise  go  to  step  (c). 

The  efficiency  of  the  algorithm  is  improved  by  scaling 
the  step  (n  +  l)’s  target  As  length  by  the  ration  of  a 
user-selected  desired  number  of  iterations  to  the  number  of 
iterations  needed  to  converge  to  step  n.  Thus,  in 
near-linear  parts  of  the  load-displacement  curve,  wide 
spacing  of  solution  points  is  allowed;  when  the  curve  rounds 
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corners  (around  limit  points)  the  method  has  to  iterate  more 
to  converge  and  hence  &s  is  reduced  until  the  curve 
straightens  out  again. 

The  inability  to  solve  for  equilibrium  at  a  limit  load 
point  is  circumvented  due  to  the  nature  of  the  stepping 
method.  The  technique  shoots  a  tangent  from  the  current 
equilibrium  point,  then  searches  an  arc  about  the  tip  for 
the  next  solution  point,  so  the  exact  limit  load  point  is 
almost  always  skipped  over.  Auxiliary  equations  can  be 
programmed  to  enable  determination  of  the  critical  point 


[49] 
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3 .  EXPERIMENTAL  METHODS 


3 . 1  Manuf actur i ng 

A  total  of  12  panels  were  fabricated  for  the  axial 
compression  experiments.  Table  3.1  breaks  down  the 
various  geometries  and  cutout  placements  used  for  the 
experimental  portion  of  the  study.  Two  additional  panels 
were  fabricated  and  used  in  the  residual  stress  test, 
described  in  more  detail  below.  All  14  panels  were  newly 
manufactured.  The  manufacturing  process  and  materials  for 
the  new  panels  were  the  same  used  by  Horban  [52] ,  Tisler 
[35] ,  and  Wilder  [53]  for  their  studies.  All  panels  were 
C-scanned  after  being  manufactured,  thus  ensuring  that 
neither  delaminations  nor  internal  defects  were  present. 

The  thickness  was  measured  at  twelve  locations  on  each 
panel,  and  the  average  thickness  for  the  panels  was  taken. 
Thickness  variation  within  each  panel  was  small  (w  6%),  and 
within  the  total  group  it  was  even  smaller  (w  3%). 

A  hole  cutting  process  was  developed  during  Tisler ’s 
work  [35]  to  help  remove  damage  caused  by  flattening  of  the 
panel  during  cutting.  A  wooden  base,  also  with  a  radiu'.  of 
12  inches,  held  the  panel  nearly  in  its  original  shape 
during  cutting.  Fiberglass  and  rubber  sheets  help  reduce 
vibration  damage,  and  a  steel  tempi  te  guided  the  router 
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Table  3.1  Experimental  Test  Plan 


Unsupported 
Vertical  Edges 

Simply  Supported 
Vertical  Edges 

8"  X  11“ 

12"  X  11“ 

8“  X  12“ 

12“  X  12" 

No  Cutout 

X 

X 

X 

X 

4“  Cutout 

X 

X 

X 

X 

4"  Cutout 

1 "  Offset 

X 

X 

X 

X 

used  for  cutting.  A  schematic  of  the  cutting  setup  is  given 
in  Fig  3.1,  along  with  photographs  in  Figs  3.2  and  3.3. 

Steel  Templale 


Fig  3.1 


Schematic  of  Cutting  Support 
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Fig  3.3 


Picture  of  Finished  Product 
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The  corners  of  the  square  cutouts  were  slightly  rounded 
(w  1/16  inch  radius).  This  was  caused  by  the  router  used  to 
cut  the  panels,  and  had  the  beneficial  effect  of  removing 
stress  singularities  at  the  corners. 

All  the  edges  of  the  panel  were  cut  the  same  way  as  was 
previously  done  [35,52,53].  NOTE:  for  those  edges  to  be 
clamped  or  simply  supported,  the  dimensions  in  the  axial  and 
circumferential  directions  were  cut  slightly  larger  (by  1/2 
inch)  to  create  a  holding  tab  for  the  experimental  fixture’s 
supports.  Uith  these  in  place,  the  effective  dimensions  in 
the  clamped  or  supported  directions  were  shorter  by  one 
inch.  It  is  crucial  that  the  upper  and  lower  edges  be 
extremely  close  to  being  parallel  to  assure  even  loading. 

The  actual  edge  cutting  procedure  is  explained  in  more 
detail  by  Horban  [52]  . 

3.2  Axial  Compression 

The  axial  compression  experimental  setup  (Fig  3.4)  and 
the  panel  clamping  devices  (Fig  3.5)  were  essentially  the 
same  as  used  by  Horban  [52] ,  Janisse  [21] ,  Tisler  [35] ,  and 
Wilder  [53]  for  their  studies  of  cutouts  and  delamination 
effects.  Initial  problems  with  the  vertical  boundary 
conditions  and  measurement  of  the  top  edge  displacement  [21] 
have  been  corrected,  as  described  by  Horban.  All  of  the 
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edges  were  held  by  specialized  fixtures,  also  explained  in 
detail  by  Horban.  The  top  and  bottom  were  clamped  (u  =  v  = 
w=R  =R  =R  =0),  while  the  edges  were  assumed  to  be 

X  y  2 

simply  supported  (u  =  R  =  free,  v=w=R=R=0)or 

X  y  z 

completely  unsupported  (u=v=w=R  =R  =R  =  free), 

X  y  z 

The  assumption  of  v  being  fixed  was  later  questioned  after 
reviewing  analytical  vs  experimental  results  (this  will  be 
discussed  in  greater  detail  in  Section  5). 

A  vertically  mounted  Linear  Variable  Displacement 
Transducer  (LVDT),  in  front  of  the  panel  in  Fig  3.4, 
measured  the  top  edge  displacement,  while  a  bank  of  LVDT’s 
behind  the  panel  measured  out  of  plane  displacement  at  15 
specified  points.  Figures  3.6  through  3.11  demonstrate  the 
LVDT  positions  for  the  8"  and  12"  panels  with  varying  cutout 
positions.  All  the  measured  radial  displacements  were  near 
the  cutout.  These  were  expected  to  be  the  largest 
displacements,  and  therefore  would  be  the  easiest  to  measure 
accurately . 

The  uniform  displacement  is  applied  from  above  with  a 
30K  lb  hydraulic  compression  machine.  Although  the  bottom 
plateen  is  moving  upward,  the  positioning  of  the  vertical 
supports  makes  t.ie  top  edge  the  true  applied  displacement. 

A  load  cell  measures  the  total  applied  load.  Uniform 
displacement  is  slowly  applied  (.05"/min)  until  a  maximum 
collapse  load  is  reached  at  which  point:  (1)  for  the  simply 
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supported  test  specimens,  the  panel  snaps  suddenly  into  a 
new  displacement  shape  and  the  load  is  reduced,  or  (2)  for 
the  unsupported  test  specimens,  the  panel  continues  to 
deform  smoothly  in  the  same  displacement  shape  it  started 
with,  while  the  load  peaks  and  then  begins  to  drop  of  slowly. 
Two  strain  gages  were  positioned  near  the  top  (Fig  3.12)  to 
ensure  the  applied  load  was  evenly  distributed.  It  is 
important  that  the  load  be  applied  evenly  by  centering  the 
force  in  the  middle  of  the  panel  to  allow  for  symmetric  load 
distribution.  For  those  panels  not  containing  cutouts,  two 
additional  strain  gages  were  positioned  at  the  geometric 
center  of  the  panel  (front  and  back),  see  Fig  3.13  for 
detail.  These  strain  gages  allowed  for  confirmation  of  the 
collapse  load  because  the  load  vs  strain  curves  diverged 
dramatically  at  collapse  (an  example  is  shown  in  Fig  3.14). 
All  the  data  from  the  strain  gages  and  LVDT’s  is  saved  on  a 
VAX  11/780  computer,  allowing  for  experimental  graphs  of  the 
load  displacement  curves,  both  for  the  top  edge  and  radial 
displacements.  For  a  detailed  discussion  of  the  axial 
compression  experimental  procedure,  including  photographs 
and  drawings  of  the  boundary  supports,  see  [21,53] . 

The  imperfection  amplitude  was  measured  using  an 
aluminum  guide  and  a  displacement  gage  [35] .  The  average 
amplitude  of  imperfections  was  found  to  be  about  the  same  as 
in  the  past  (w  0.0060  inches).  Note  that  this  is 
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Strain  Gages 
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NOTE:  FOR  SOLID  PANELS,  2  STRAIN 
GAGES  ARE  PLACED  AT  CENTER 
OF  PANEL  ON  BOTH  SIDES. 


Fig  3.13  Strain  Gage  Positions,  Panels  Without  Cutouts 
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approximately  the  thickness  of  one  ply  and  is  substantially 
less  than  the  panel  thickness  of  0.042  inches. 

3.3  Residual  Stress  Test 

For  the  determination  of  residual  stress  created  by 
cutout  removal,  16  strain  rosettes  are  mounted  along  the 
edge  of  the  cutout  location  (Figs  3.15  and  3.16).  These 
rosettes  are  mounted  prior  to  removing  the  cutout.  Because 
of  the  type  of  router  used,  the  rosettes  are  placed  within 
0.5  inches  of  the  cutout  edge  only.  Strain  in  the  panel 
while  free  from  any  load,  is  recorded  to  create  a  baseline 
of  strain  data  which  determines  the  amount  of  residual 
stress  created  by  removing  the  cutout  (Fig  3.17).  The  panel 
is  placed  in  the  cutting  setup,  and  again  strain  is  recorded 
generating  a  baseline  of  strain  data  while  the  panel  is 
strapped  into  the  cutting  setup.  After  the  router  punched  a 
hole  into  the  panel,  strain  is  recorded  (Fig  3.18),  and 
strain  is  recorded  after  each  side  of  the  cutout  template  is 
completed  by  the  router .  Once  the  cutout  is  completely 
removed,  strain  is  recorded  again  before  the  panel  is 
removed  from  the  cutout  setup  (Fig  3.19).  The  panel  is 
removed  from  the  cutout  assembly  and  placed  such  that  no 
loading  is  applied.  Then  strain  is  measured  for  a  final 
time,  for  a  total  of  nine  strain  measurements.  Residual 
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stress  is  determined  from  any  changes  in  the  measured  strain 
from  the  first  and  last  strain  measurements.  The  other 
measurements  are  recorded  as  an  attempt  to  give  a  complete 
history  of  the  cutout  process  and 


Fig 
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Strain  Rosette  Positions,  Front  View 
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Fig  3.16 


Strain  Rosette  Positions,  Rear  View 
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aid  in  determining  any  anomalies  that  may  be  created  by 
removing  the  cutout.  The  change  in  strain  is  attributed  to 
cutting  out  a  portion  of  the  curved  composite  panel  .  A 
better  approximation  of  the  residual  stress  distribution  is 
generated  than  in  [35]  but  it  by  no  means  an  exact 
distribution . 

When  a  curved  panel  is  cured  from  an  elevated 
temperature,  a  residual  stress  forms.  These  stresses  are 
considered  to  have  two  parts.  The  first  is  an  in-plane 
stress,  with  no  warping  involved  if  the  panel  is  symmetric 
and  balanced.  This  stress  can  be  easily  calculated  by 
assuming  a  flat  plate  with  the  ply  lay-up  used  for  this  work. 
The  second  part  comes  from  the  curvature.  As  the  curved 
panel  is  cured,  the  uneven  lengths  of  the  inner  and  outer 
plies  create  a  bending  stress.  This  is  referred  to  as  the 
bending  residual  stress. 

The  stress  caused  by  curing  a  flat  composite  plate  with 
the  ply  lay-up  [0 ,-45 ,+45 ,90]  is  estimated  using  a  curing 
temperature  change  of  300 °F.  This  flat  plate  stress  is 
found  to  be  %  5  ksi .  A  bending  stress  also  comes  about 
as  a  result  of  the  panel’s  curvature,  as  discussed  above. 

If  a  change  occurs  within  the  flat  plate  residual  stress, 
the  change  should  be  the  same  in  all  directions  since  the 
panel  is  quasi-isotropic.  If  a  change  takes  place  in  the 
bending  residual  stress,  the  change  will  predominantly  be  in 
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the  circumferential  direction,  with  limited  change  in  the 
axial  direction.  By  comparing  the  change  in  strain  in  the 
circumferential  and  axial  directions,  one  can  tell  which 
component  of  the  residual  stress  is  affected  by  the  cutout. 
Finally,  changes  in  stress  due  to  changes  in  strain  are 
compared  to  the  stresses  present  during  axial  compression  to 
see  if  residual  stress  change  warrents  further  study. 


FINITE  ELEMENT  MODELING 


4.1  Panel  Properties  and  Analytical  Boundaries 


Once  placed  within  the  experimental  setup,  all  panels 
have  a  radius  of  12“  but  axial  and  circumferential 
dimensions  varied.  These  are  discussed  in  Table  3.1  of 
Section  3.  Essentially,  there  are  two  classes  of  panels. 

One  class  is  clamped  on  the  top  and  bottom  and  had  free 
vertical  edges.  Of  this  class,  panels  is  either  8“  or  12“ 
in  the  circumferential  direction  and  11"  in  the  axial 
direction.  Each  circumferential  dimension  (8"  or  12")  has 
three  panels  within  its  sub-class:  one  panel  has  no  cutout, 
one  panel  has  a  4"  cutout  centered  within  the  panel,  and  the 
final  panel  has  a  4“  cutout  offset  one  inch  to  measure 
eccentricity  effects  of  the  cutout.  The  other  class  is 
clamped  on  the  top  and  bottom  and  has  simply  supported 
vertical  edges.  This  class  is  broken  down  in  the  exact 
manner  as  the  first  class.  The  panel  material  properties, 
dimensions,  and  sign  conventions  in  Fig  4.1.  These  material 
properties  are  determined  experimentally  through  tensile 
test  coupons  and  are  used  for  the  finite  element  model. 

For  the  finite  element  models,  as  mentioned  before,  two 
boundary  conditions  are  compared.  In  both  models,  the 
bottom  boundary  is  clamped  (u=v=w=R=R=R=0). 


Material:  Hercules  AS4  -  3501  Graphite  Epoxy 

Radius,  R:  12irx:hes 

Length,  L 1 1  or  12  inches 

Ciraimterentlal  Length:  8^  or  1 2“  ®  oo 

Ply  Layup:  [0,  -  45,45,90^ ,  angles  as  indicated 

Elastic  Moduli:  E^=  20.461  xio’psi 

E2=  1.3404x10*  psi 

Shear  Modulus:  0^2*  -^638  x  10*  psi 
Poisson’s  Ratioe:\^j- .301  .0206 

Fig  4.1  Properties  for  Composite  Panels 
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And  in  both  models,  the  top  has  a  degree  of  freedom,  u, 
free.  For  the  one  class,  the  vertical  edge  boundaries  are 


free  (u=v=w=R  =R  =R  /O).  For  the  other  class, 

U  V  w 

the  vertical  edge  boundaries  are  assumed  simply  supported, 
both  vertical  rotation  and  movement  (R  and  u)  are  assumed 

U 

free . 


Both  Janisse  [21]  and  Tisler  [35]  compared  simply 
supported  panels  and  varied  whether  the  degree  of  freedom, 

V,  is  fixed  or  free.  Janisse  [21]  found  this  variation  to 
have  an  almost  negligible  effect  for  a  panel  with  a  2" 
cutout  (3%  difference),  but  anticipated  the  difference  to 
be  more  pronounced  for  larger  cutouts.  Because  the 
difference  is  so  small  for  2"  cutouts,  he  did  not  have  to 
consider  whether  the  experimental  edge  boundary  conditions 
were  v  free  or  v  fixed.  Tisler  [35]  on  the  other  hand  did 
consider  this  variation  for  a  panel  with  a  4"  cutout  (18% 
difference)  and  found  it  to  have  a  significant  effect  on  the 
panel’s  global  collapse  load.  Because  of  some  disparity 
between  early  analytical  (STAGS)  and  experimental  results, 
those  models  not  containing  cutouts  are  run  with  the  degree 
of  freedom,  v,  free,  and  fixed,  to  determine  which  boundary 
condition  better  approximated  the  true  experimental 
conditions . 
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4.2  Grid  Selection 


Based  on  thorough  research  into  the  previous  works  of 
Lee  [22],  Tisler  [35],  and  Egan  [37],  the  411  element  (Fig 
4.2)  proved  to  be  the  best  choice  for  use  in  the  nonlinear 
large  displacement/moderate  rotation  analysis  in  STAGS. 
Unfortunately,  the  lack  of  a  transition  element  made  grid 
selection  important.  Again,  careful  analysis  of  Lee’s  and 
Tisler ’s  work  determined  a  refined  grid  utilizing  a  0.5  inch 
element  around  the  cutout  provided  the  best  results.  Hence 
grids  similar  to  those  shown  in  Figs  4.3  through  4.6  are 
used.  These  grids  are  slightly  more  refined  than  Tisler ’s. 
This  is  due  to  the  large  displacements  (w  .02“)  and 
rotations  (w  30°)  experimentally  measured,  and  analytically 
determined,  for  the  panels  with  unsupported  vertical  edges. 

It  must  be  noted  that  these  refined  grid  selections  require 
large  amounts  of  computer  time  on  the  VAX  computer  systems 
(up  to  90000  CPU  seconds)  in  order  to  reach  collapse  load. 

Based  on  Dennis’  work  [41]  ,  the  36  degree  of  freedom 
(DOF)  element  (Fig  4.7)  proved  to  be  the  best  choice  for  use 
in  the  nonlinear  analysis  in  SHELL.  Again  grid  selection  was 
important,  but  another  factor  presented  itself.  In  order  to 
minimize  the  amount  of  CPU  run-time  needed  to  converge  to  an 
accurate  solution  ( for  the  number  of  models  that  are 
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considered),  it  was  decided  to  use  only  half  size  models 
(examples  are  shown  in  Figs  4.8-10)  on  the  computer 
mainframes  available  (Elexi  6400  and  VAX  8550)  to  converge 
to  a  solution  using  SHELL.  Due  to  the  asymmetrical 
displacement  pattern  that  exists  at  the  point  of  collapse, 
it  is  recognized  these  models  will  generate  solutions  with 
larger  stiffnesses  (thereby  requiring  larger  collapse 
loads).  An  attempt  was  made  to  determine  the  percentage 
increase  over  the  full  model  by  using  a  comparison  between 
full  and  half  models  using  a  coarser  grid  (see  Figs  4.11  and 
4.12).  It  is  recognized  that  this  comparison  is  in  no  way 
an  accurate  method  in  determining  how  much  increase  will 
exist  for  each  model.  This  is  just  an  attempt  in  arriving 
to  a  rough  estimate  as  to  a  range  of  increase  that  could  be 
applied  to  the  half  models  and  thereby  reduce  the  artificial 
stiffness  created  by  the  model(s).  The  artificial  stiffness 
arises  due  to  the  fact  that  the  displacement  pattern  is 
asymmetrical.  This  asymmetry  is  created  by  the  and  the 
D  bending  coefficients  (which  in  turns  affects  the  M 

26  xy 

moment)  coming  in  to  play  at  the  collapse  of  the  panel. 

Once  the  percentage  increase  is  determined,  it  is  applied  to 
the  half  model  solutions  as  a  means  of  removing  the 
artificially  induced  stiffness. 


STAGS  Grid  for  8"  Panel,  Centered  Cutout 
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1/2" 

VTU 


Fig  4.6  STAGS  Grid  for  12"  Panel,  Centered  Cutout 
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5.  RESULTS  AND  DISCUSSION 


The  main  study  area  of  this  thesis  is  the  correlation 
of  analytical  results  with  experimentally  derived  results. 
This  correlation  considers  varying  panel  geometry,  boundary 
conditions,  geometric  imperfection  (cutout)  effects 
including  eccentricity  effects,  and  through-the-thickness 
shear  effects.  Thus,  in  order  to  establish  the  critical 
characteristics  of  the  effect  that  a  true  geometric 
imperfection  (cutout)  has  on  a  composite  panel,  a  study  was 
carried  out  incorporating  the  previously  mentioned 
comparisons.  The  results  of  these  comparisons  are  discussed 
below.  The  discussion  includes  the  analytical  results  from 
STAGS  and  SHELL,  how  these  results  correlate  with 
experimentally  derived  results,  and  a  brief  discussion  of 
determining  residual  stress  experimentally. 

5 . 1  Analytical  Results  -  STAGS 

Before  beginning  this  discussion,  it  should  be 
reiterated  what  conditions  are  being  started  with,  since 
several  parameters  were  varied.  One  begins  with  the 
assumption  of  a  12"  wide  panel,  no  cutout,  no  initial  radial 
imperfections  included,  clamped  along  the  top  and  bottom, 
and  simply  supported  along  the  vertical  edges  (assuming  v 
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is  fixed  along  the  vertical  edges).  Then  panel  width, 
boundary  conditions,  geometric  imperfections  (cutouts),  and 
eccentricity  of  cutout  location  are  mentioned  at  the 
start  of  each  model's  discussion. 

Based  on  previous  work  [35],  it  is  assumed  the  loading 
condition  is  a  top  edge  uniform  displacement.  During  the 
experimental  collapse  tests,  a  top  edge  uniform  displacement 
was  applied,  not  a  uniform  load,  to  the  panels  through  the 
test  fixture.  Tisler  [35]  determined  that  a  uniform 
displacement  load  condition  for  STAGS  best  matches 
experimental  results.  With  a  uniform  displacement,  applied 
load  is  measured  by  summing  the  total  axial  force  along  the 
top.  This  is  done  by  summing  the  equilibrium  load  values 
for  each  node  along  the  top  of  the  panel  (Fig  5.1).  This 
load  is  far  from  being  a  constant  distribution  when  a  large 
cutout  (11-16%  reduction  in  surface  area)  is  present  as  seen 
in  Fig  5.2.  By  summing  the  total  axial  force  analytically, 
one  can  make  comparisons  with  the  total  load  measured 
experimentally. 

5.1.1  V  Free  vs  V  Fixed 

During  current  and  past  cutout  experiments  [24,35],  it 
became  apparent  that  the  assumed  vertical  boundary  condition 
of  V  free  was  not  entirely  consistent.  Teflon  tape 
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4  or  6  inches 


Pig  5.1  Equilibrivun  Load  Distribution  along  Top  of  Solid 
Panel  at  Collapse 


□ 

□ 

□ 

□ 

□ 

□ 

□  □ 

Fig  5.2  Equilibrium  Load  Distribution  along  Top  of  4 
Cutout  Panel  at  Collapse 
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protected  the  edges  and  corners  of  the  panel  from  being 
damaged  by  the  experimental  test  setup.  Although  there  was 
scarring  of  this  tape  during  the  tests  (indicating  the 
circumferential  (v)  displacement  was  somewhat  restricted), 
it  was  not  clear  if  indicated  v  was  completely  fixed  along 
the  vertical  edges.  A  comparison  was  carried  out  for  the 
same  model  with  v  fixed  and  v  free. 

The  added  restriction  of  v  fixed  along  the  vertical 
edge  allows  for  more  of  the  load  to  be  carried  along  the 
vertical  edges,  so  one  expects  the  panel  with  v  fixed 
to  carry  more  load  before  it  collapsed.  The 
load-displacement  curves  for  a  solid  panel  are  shown  in  Figs 
5.3  and  5.4,  comparing  v  free  to  v  fixed.  It  is  seen,  as 
expected,  the  fixed  boundary  condition  results  in  a  higher 
collapse  load  (about  10%)  and  top  edge  displacement.  This 
collapse  load  increase  is  much  greater  than  the  3%  increase 
noted  by  Janisse  [21]  for  2”  cutouts  but  less  than  the  18% 
increased  noted  by  Tisler  [35]  for  5"  cutouts.  As  Tisler 
noted,  the  radial  displacements  for  the  v  fixed  are  greater 
than  for  v  free  (Pigs  5.5-8).  With  v  fixed,  most  of  the 
bending  and  load  carrying  occurs  at  the  vertical  edges  of 
the  panel,  while  they  are  more  evenly  distributed  for  the  v 
free  condition. 

When  compared  to  experimental  results  for  collapse  load 
and  displacement  of  the  12"  panel  with  no  cutout  (Fig  5.3) 
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Total  Applied  Load  (lbs) 


Fig  5.3 


Load-Displacement,  v  Free  and  Fixed,  12"  Panel 
Experimental  Results  Compared 
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Total  Applied  Load  (lbs) 


Fig  5.4  Load-Displacement,  v  Free  and  Fixed, 
Experimental  Results  Compared 


V  Fixed 

V  Free 
Experimental 


8"  Panel 
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Fig  5.5  Radial  Displacement  Contour,  12”  Wide  Panel,  No 
Cutout,  V  Free,  Total  Load  =  4754  lbs  (100%  of 
Collapse  Load). 


Fig  5.7  Radial  Displacement  Contour,  8"  Wide  Panel,  No 
Cutout,  V  Free,  Total  Load  =  3174  lbs  (100%  o£ 
Collapse  Load) . 


and  the  8"  panel  with  no  cutout  (Fig  5.4),  it  appeared  that 
the  two  conditions  form  an  upper  and  lower  bound  for  the 
experimental  results.  From  a  physical  viewpoint,  based  upon 
scarring  of  the  teflon  tape,  it  seemed  better  to  use  the  v 
fixed  condition  to  model  the  response.  Although  the  author 
believes  that  the  actual  experimental  boundary  conditions 
probably  lies  somewhere  between  the  two  extremes  of  fixed 
and  free,  all  further  analytical  comparisons  are  made 
assuming  that  v  is  fixed  along  the  vertical  boundaries  for 
the  simply  supported  models.  Because  of  this  simplifying 
assximption,  the  STAOS  collapse  loads  are  expected  to  be 
slightly  higher  than  the  experimental  results. 

5.1.2  Surface  Imperfections 

Janisse  showed  the  greatest  reduction  in  collapse  load 
for  uncut  panels  by  modeling  imperfections  as  three  half 
sine  waves  in  both  the  axial  and  circumferential  directions 
[21].  It  is  felt  that  this  imperfection  modeling  more 
closely  approximates  the  experimental  results  (for  uncut 
panels  only)  in  part  because  the  radial  displacement  pattern 
at  collapse  somewhat  approximates  this  (refer  Figs  5.5-8). 

It  is  noted  that  v  is  now  assvuned  fixed  along  the  vertical 
boundaries . 

Imperfection  models  using  one,  two,  or  four  half  sine 
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waves  all  yielded  higher  collapse  loads  (from  as  little  as 
15%  higher  up  to  75%  higher)  than  the  perfect  panel  model . 
The  hypothesis  is  that  since  the  buckled  shape  somewhat 
approximates  three  half  sine  waves  in  both  the  radial  and 
axial  directions,  it  requires  less  energy  to  move  through 
an  imperfection  of  a  similar  shape  to  the  collapsed 
shape  of  the  panel,  than  for  the  other  imperfection  models 
used.  Hence,  any  imperfection  models  included  in  the 
following  discussion  consisted  of  the  three  half  sine  waves 
with  an  maximum  amplitude  of  0.005"  (w  one  ply  thickness). 

With  the  imperfections  included,  the  computer  time 
required  to  reach  collapse  load  was  at  least  doubled  (from 
35000  CPU  seconds  to  over  72000).  The  final  reduction  in 
load  turned  out  to  be  small  (Table  5.1,  Figs  5.9-12), 
between  1-3%.  This  is  especially  true  when  considered 
with  the  very  large  reductions  seen  by  Janisse  [24]  for 
solid  panels  (no  cutout).  For  an  equivalent  ply  lay-up  and 
imperfection,  his  uncut  panel  with  the  same  imperfection 
pattern  collapsed  at  less  than  25%  of  the  bifurcation 
collapse  load.  Note:  for  a  panel  with  no  cutouts  or 
imperfections,  the  bifurcation  solution  is  essentially  the 
same  as  the  nonlinear  solution  {*»  3%  variance).  By 
modeling  more  precisely  the  imperfections  typically  seen  in 
panels  without  cutouts,  Starnes  [30]  saw  collapse  load 
reductions  of  only  about  10%.  Since  collapse  load  of  panel 


with  cutouts  was  much  less  sensitive  to  imperfections  in  the 
worst  case  form,  it  seems  reasonable  to  assume  sensitivity 
will  also  be  less  for  other  imperfection 
patterns  such  as  those  mentioned  above. 

When  an  imperfection  model  was  tried  with  a  cylindrical 
panel  containing  a  large  cutout  it  was  found  to  be  less 
sensitive  to  imperfections  than  the  same  panel  with  no 
cutout.  Previous  work  by  Tisler  [35]  seemed  to  bear  this 
out.  However,  this  was  not  the  primary  thrust  of  his  work 
and  therefore  does  not  give  much  depth  to  this  particular 
analysis.  Unfortunately,  no  other  work  on  cylinders  with 
both  large  cutouts  and  imperfections  was  found  to  verify 
these  results. 


Table  5.1  Results  by  Comparing  Perfect  vs  Imperfect  Panels 


Collapse  Load,  lbs 


Panels 

Perfect 

Imperfect 

Reduction 

12" 

SS 

5060 

5940 

2.2% 

8" 

SS 

3761 

3700 

1.6% 

12" 

US 

2298 

2260 

1.7% 

8" 

US 

1382 

1350 

2.3% 

MOTE:  SS  =  Simply  Supported,  US  =  Unsupported 


Total  Applied  Load  (lbs) 


Fig  5.9  Load-Displacement,  12"  Simply  Supported 
Panel,  Comparing  Perfect  vs  Imperfection 
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Fig  5.10  Load-Displacem 
Panel ,  Compari 


Total  App 


Fig  5.11  Load-Displacement,  12"  Unsupported  Panel, 
Comparing  Perfect  vs  Imperfection 


Total  Applied  Load  (lbs) 


Perfect  Panel 
Imperfection 


Fig  5.12 


Load-Displacement,  8"  Unsupported  Panel, 
Comparing  Perfect  vs  Imperfection 


5- 


It  should  be  noted  at  this  point  that  for  the  solid 
panels  with  unsupported  edges  had  extreme  problems 
converging  to  a  solution  with  the  STAGS  QUAF  411  element. 
However,  when  the  QUAF  410  element  was  tried,  along  with  the 
imperfections,  STAGS  did  converge  to  solutions  that  were 
reasonable  when  compared  to  the  experimental  results.  The 
author  will  discuss  this  in  greater  detail  in  Section  5.3. 

The  STAGS  radial  displacement  collapse  patterns  (Figs 
5.13-16)  and  isoparametric  views  (Figs  5.9-12)  were  much 
different  for  panels  with  imperfections  than  for  panels 
without  at  the  collapse  load.  This  imperfection  pattern 
alluded  to  earlier  in  the  discussion  is  not  the  true 
experimental  imperfection  pattern,  and  the  resulting 
collapse  shapes  are  much  different.  For  those  models  with 
an  surface  imperfection  model,  a  trough  (aligned  with  the 
-45**  outer  angle  ply)  is  distinctly  formed  whereas  in  the 
perfect  model,  the  trough  is  barely  distinguishable  and  is 
nearly  perpendicular  (aligned  with  the  0®  ply).  The 
experimental  radial  displacements  were  very  similar  to  that 
predicted  for  panels  with  no  imperfections.  This  indicates 
that  initial  imperfections  do  not  control  the  displacement 
pattern  to  the  extent  seen  in  the  analytical  imperfection 
case.  However,  as  evidenced  in  Figs  5.9-12,  the  initial 
imperfections  do  reduce  the  panel's  stiffness  (hence  its 
ability  to  withstand  additional  energy  created  by  increased 
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Fig  5.13  Radial  Displacement  Contour,  12"  Wide  Panel, 
No  Cutout,  Simply  Supported,  No  Imperfection, 
Total  Load  =  5060  lbs  (100%  of  Collapse  Load) 


Fig  5.14  Radial  Displacement  Contour,  12"  Wide  Panel, 
No  Cutout,  Simply  Supported,  Imperfection, 
Total  Load  =  4899  lbs  (100%  of  Collapse  Load) 
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Fig  5.15  Radial  Displacement  Contour,  8"  Wide  Panel, 

No  Cutout,  Simply  Supported,  No  Imperfection, 
Total  Load  =  3761  lbs  (100%  of  Collapse  Load) 
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Pig  5.16  Radial  Displacement  Contour,  8”  Wide  Panel, 

No  Cutout,  Simply  Supported,  Imperfection, 
Total  Load  =  3749  lbs  (100%  of  Collapse  Load) 


loading)  over  the  duration  of  the  load-displacement  curve. 

It  is  important  to  notice  that  the  unsupported  panels  are 
significantly  affected  by  the  initial  imperfections. 

Although  the  analytical  imperfection  model  used  does  not 
accurately  reflect  the  experimental  imperfections,  it  does 
demonstrate  the  inherent  reduction  in  stiffness  exhibited  in 
the  panel  experimentally.  Also,  in  the  case  of  the 
unsupported  panels,  the  imperfection  model  does  indicate 
increased  bending  in  the  panel's  response  once  it  has 
entered  the  nonlinear  region  of  loading.  The  reduction  in 
collapse  load  due  to  initial  imperfections  is  present  for 
all  the  panels  shown  in  the  above  figures.  In  the  case  of 
the  simply  supported  panels,  the  analytical  imperfection 
model  demonstrated  a  small  change  in  collapse  load  and 
stiffness  (w  5-7%).  In  the  unsupported  panels  the 
analytical  imperfection  model  showed  a  small  change  in 
collapse  load  (m  5%),  but  a  marked  decrease  in  stiffness 
with  a  significant  increase  in  bending  (as  much  as  a  35% 
increase  in  M  ) . 

xy 

Because  of  the  small  change  in  collapse  load  (less  than 
5%)  at  the  cost  of  a  substantial  increase  in  computer  time 
(over  75000  CPU  seconds)  and  poor  radial  displacement 
predictions  (it  was  not  practical  to  use  the  methods  needed 
to  get  the  more  accurate  initial  imperfection  shape  [30]), 
it  was  decided  to  perform  the  analytical  checks  without 
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imperfections  for  the  panels  with  cutouts.  It  is  assumed 
based  on  previous  work  [35]  and  the  few  computer  models 
attempted  in  this  effort  that  the  4”  cutout  panels  are 
relatively  insensitive  to  initial  imperfections.  The  author 
feels  that  the  large  geometric  imperfections  for  these 
panels  dominate  the  panels  response  over  the  initial 
imperfections.  As  in  assuming  a  fixed  circumferential 
boundary  condition  (v),  the  STAGS  solution  that  assumes 
no  initial  imperfections  should  estimate  higher  collapse 
loads  than  is  seen  experimentally.  Only  the  panels  with  no 
cutouts  contain  imperfection  models,  primarily  due  to  the 
relatively  small  cost  of  running  these  models  (only  an 
increase  of  about  10000-15000  CPU  seconds)  and  the  marked 
change  in  stiffness  and  bending  due  to  the  presence  of  the 
initial  imperfections. 

5.1.3  Geometric  Imperfections  (Cutouts)  &  Eccentricity 
Effects 


One  of  the  major  thrusts  of  this  effort  is  studying  the 
collapse  characteristics  of  cylindrical  composite  shells  due 
to  the  effects  of  large  geometric  imperfections  (cutouts) 
and  the  effect  of  cutout  eccentricity. 

Load-displacement  curves  were  generated  for  all  twelve 
of  the  models  mentioned  in  Section  3.1,  Table  3.1  and  are 
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shown  in  Figs  5.17-20  in  composite  form  (meaning  curves  are 
shown  together  for  the  same  boundary  conditions  and  panel 
width).  It  should  be  noted  that  v  is  assumed  fixed  for  all 
simply  supported  panels  and  that  all  panels  with  no  cutouts 
do  contain  an  initial  imperfection  model  of  three  half  sine 
wave  in  both  the  radial  and  circumferential  directions.  The 
maximum  amplitude  of  these  sine  waves  is  one  ply  thickness 
(0.005") . 

By  comparing  values  between  the  same  class  of  panels, 
it  is  clear  the  effect  cutouts  have  on  the  panel  collapse 
characteristics.  It  is  also  evident  the  effect  eccentricity 
has  on  the  panel (s).  In  the  case  of  the  12",  simply 
supported  panels,  the  4"  cutout  decreases  the  collapse  load 
by  51%.  However  cutout  eccentricity  actually  increases  the 
collapse  load  by  a  small  amount  (3%)  over  the  centered 
cutout.  Both  panels  with  cutouts  show  an  extensive  amount 
of  nonlinearity  when  compared  with  the  panel  that  has  no 
cutout.  As  Tisler  [34],  Hermsen  [23],  and  Janisse  [21] 
noted,  this  nonlinearity  is  caused  by  the  large  stress 
(axial  and  bending)  redistribution  occurring  with  large 
cutouts.  The  difference  between  the  nonlinear  collapse 
loads  and  buckling  loads  obtained  by  simpler  and  more 
efficient  bifurcation  analysis  is  seen  to  change  with  cutout 
size  and  location  (Fig  5.21).  As  expected,  for  no  cutouts 
with  and  no  initial  imperfections,  the  nonlinear  and 
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No  Cutout  flmperfect' 
Cutout,  Centered 
Cutout,  Offset 


Pig  S.17  Load-Displacement,  12”  Simply  Supported 
Panel,  Comparing  Geometric  Imperfections 


Fig  5.18  Load-Displacement,  8"  Simply  Supported 

Panel,  Comparing  Geometric  Imperfections 
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Fig  5.19  Load-Displacement,  12"  Unsupported 

Panel,  Comparing  Geometric  Imperfections 


# 


bifurcation  values  are  the  same.  The  2"  cutout  does  not 
have  much  stress  redistribution,  and  the  bifurcation  value 
is  higher  than  the  nonlinear  collapse  load  [21].  The 
opposite  occurs  for  the  4"  cutouts,  with  the  bifurcation 
values  becoming  lower  than  the  nonlinear  results.  The  lower 
bifurcation  values  for  the  4”  cutouts  is  expected,  since 
linear  bifurcation  methods  do  not  consider  the 
redistribution  of  moments  and  stresses.  Instead,  the 
bifurcation  method  assumes  that  the  stresses  and  moments 
increase  linearly  from  a  small  initial  load.  Fig  5.22  shows 
how  the  vertical  stress  profile  (N  )  in  the  panel  with  the 

X 

4"  cutout,  centered,  changes.  This  N  profile  is  taken 

X 

along  the  left  side  of  the  cutout.  The  stress  is  greatest 
near  the  cutout  at  the  beginning  of  loading,  and  shifts 
toward  the  vertical  edges  while  loading  is  increased.  At 
collapse,  most  of  the  load  is  carried  along  the  vertical 
edges  during  loading.  This  shifting  of  the  load  to  the 
supported  edges  allows  for  a  much  higher  collapse  load  than 
is  predicted  by  the  linear  bifurcation  method.  When  the 
cutout  is  moved  closer  to  one  side,  this  redistribution 
occurs  even  more  quickly.  For  the  4"  centered  cutout,  the 
redistribution  was  essentially  accomplished  at  80%  of  the 
collapse  load.  For  the  4"  cutout,  offset  1",  the 
redistribution  was  accomplished  at  approximately  70%  of  the 
collapse  load.  This  may  indicate  the  reason  why  the 
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Bifurcation/Nonlinear  Collapse  Load 


Fig  5.23  Bifurcation/Nonlinear  Collapse  Load  for 
Varying  Cutout  Sizes  &  Location 
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Nx,  Ibs/alament  width 


Fig  5.24  Values  o£  N  Along  Side  of  Cutout,  12"  Simply 
Supported  Panel 
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eccentricity  of  the  cutout  leads  to  a  slight  increase  in  the 
total  collapse  load  of  the  panel. 

The  same  trend  is  noticed  for  the  other  three  classes 
of  panels.  In  the  case  of  the  8"  simply  supported  panel, 
the  4"  centered  cutout  reduces  the  total  collapse  load  by 
43%.  As  in  the  12"  simply  supported  panels,  eccentricity  of 
the  cutout  increases  collapse  load  by  3%.  The  smaller 
reduction  in  collapse  load  is  attributed  to  the  simply 
supported  edges  being  closer  to  the  cutout,  thereby 
influencing  the  stress  redistribution.  Fig  5.23  shows  how 
the  vertical  stress  profile  for  the  8"  panel  with  a  centered 
cutout  changes.  As  expected,  the  redistribution  follows  the 
pattern  shown  by  the  12"  panel  with  the  off-center  cutout. 

As  the  simply  supported  vertical  edges  are  moved  closer  to 
the  cutout,  the  stress  redistribution  occurs  at  a  lower 
percentage  of  the  total  collapse  load.  For  the  8"  panel 
this  occurred  at  approximately  65%  of  the  collapse  load. 

And  for  the  8"  panel  with  the  offset  cutout,  the 
redistribution  occurred  at  approximately  50%  of  the  collapse 
load. 

For  the  12"  unsupported  panels,  the  presence  of  a  4" 
centered  cutout  reduced  the  total  collapse  load  by  46%.  But 
when  the  cutout  was  offset  1",  the  collapse  load  is  reduced 
only  by  40% 

The  8"  unsupported  panels  showed  extreme  sensitivity  to 


Nx,  lb$/elemeni  width 


Fig 


5.23  Values  of  Along  Side  of 
Supported  Panel 
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the  presence  of  cutouts.  The  4"  centered  cutout  reduced  the 
total  collapse  load  by  77%,  and  when  the  cutout  was  offset 
1",  the  collapse  load  is  reduced  by  70%. 

5.1.4  Boundary  Conditions 


I 


I 


Another  major  thrust  of  this  effort  is  studying  the 
effect  of  changing  the  boundary  conditions  along  the 
vertical  edges  of  the  composite  panels.  Load-displacement 
curves  were  generated  in  comparing  six  categories  of  panels 
based  on  panel  geometry  (width)  and  the  existence  of 
geometric  imperfections.  These  are  shown  in  Figs  5.24-29. 
In  every  circumstance,  the  removal  of  the  vertical  simple 
supports  resulted  in  a  significant  decrease  of  the  total 
collapse  load.  As  stated  before,  all  solid  panels  contain 
initial  imperfection  models  of  three  half  sine  waves  with 
amplitude  of  one  ply  thickness,  v  is  assumed  to  be  fixed 
for  all  simply  supported  panels.  It  is  expected,  based  on 
these  initial  assumptions,  that  STAGS  will  converge  to  a 
solution  with  more  stiffness  than  is  seen  experimentally. 

Table  5.2  shows  the  tabulated  comparison  for  the 
analytical  results  from  STAGS.  Fig  5.30  shows  the 
comparison  between  simply  supported/unsupported  panels  for 
each  category.  On  the  average,  for  the  12" 
panels,  the  reduction  is  50-53%.  For  the  8"  panels,  the 
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Fig  5.25  Load-Displacement,  12”  Panels,  4"  Cutout,  Centered 
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Fig  5.27  Load-Displacement,  8”  Panels,  No  Cutout 


5-39 


loial  Applied  Lood  \^lbs 


Fig  5.28  Load-Displacement,  8"  Panels,  4"  Cutout,  Centered 
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Pig  5.29  Load-Displacement,  8”  Panels,  4"  Cutout,  1”  Offset 
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reduction  is  substantial,  between  63-85%. 


Table  5.2  Panel  Results  by  Comparing  Boundary  Conditions 


Panels 

SS 

US 

Reduction 

12",  NC 

4950 

2260 

54% 

12",  CC 

2521 

1186 

53% 

12",  OC 

2561 

1284 

50% 

8",  NC 

3700 

2260 

63% 

8",  CC 

2178 

323 

85% 

8",  OC 

2465 

460 

81% 

NOTE:  NC  =  No  Cutout,  CC  =  Centered  Cutout 
OC  =  Off  center  Cutout,SS  =  Simple 
Support,  US  =  Unsupported 

Of  all  the  panels  considered,  the  8"  wide  panels  showed 
and  extreme  reduction  when  the  vertical  supports  are 
removed.  Combine  this  with  the  presence  of  a  geometric 
imperfection  and  the  reduction  is  almost  complete.  A  quick 
inspection  of  the  data  thus  far  indicates  the  8"  panels 
(with  a  geometric  imperfection  present)  are  behaving  as  if 
they  were  composed  of  two  distinct  columns  that  are  joined 
by  some  material  at  the  top  and  bottom.  This  in  turn, 
indicates  these  panels  are  no  longer  behaving  as  shells  but 
as  thin,  wide  plates  or  (if  the  plates  are  oriented 


Collapse  Load  (lbs) 


correctly)  as  colunms.  However,  the  author  was  unable  to 
model  this  assumption  using  some  simplified  theoretical 
relations  from  [4,44].  By  using  the  simplified  column  and 
plate  relations  (Euler  column,  Donnell  approximation)  the 
error  between  theory  and  the  analytical  results  was  an 
average  of  40%  with  one  model  approaching  70%. 

The  only  positive  aspect  of  this  assumption  is  that  it 
may  indicate  why  a  panel,  with  an  offset  cutout,  yields  a 
higher  collapse  load.  When  offsetting  a  cutout  by  1",  this 
creates  one  column  that  is  wider  than  in  the  case  of  the 
panel  with  a  centered  cutout,  and  may  play  a  role  in 
increasing  the  panel's  overall  stiffness.  However,  as  is 
stated  above,  the  panel's  response  is  too  complicated  to  be 
explained  by  relatively  simple  column  and  plate  relations. 
Unfortunately,  there  exist  no  literature  that  addresses 
axially-loaded,  cylindrical  composite  shells  with  large 
cutouts  having  unsupported  vertical  edges.  Therefore,  to 
gain  a  further  understanding  of  this  complex,  nonlinear 
problem,  further  analytical  and  experimental  research  must 
be  accomplished. 

5.1.5  Panel  Geometry 

The  final  major  portion  of  this  effort  is  considering 
boundary  effects  upon  cylindrical  composite  shells  which 


have  large  geometric  imperfections.  Load-displacement 
curves  were  generated  in  comparing  six  categories  of  panels. 
These  categories  compared  panel  width  while  keeping  the 
other  parameters  (boundary  conditions,  geometric 
imperfections)  constant.  These  are  shown  in  Figs  5.31-36.  In 
every  circumstance,  the  reduction  in  panel  width  for  the 
same  boundary  conditions  and  the  same  geometric  imperfection 
(or  lack  thereof)  yielded  a  reduction  in  total  collapse 
load.  As  stated  before,  all  panels  with  no  geometric 
imperfections  contained  an  initial  imperfection  model  of 
three  half  sine  waves  in  both  circumferential  and  radial 
directions  with  an  amplitude  of  one  ply  thickness,  v  is 
assumed  to  be  fixed  for  all  simply  supported  panels.  It  is 
expected,  based  on  these  initial  assumptions,  that  STAGS 
will  converge  to  a  solution  with  a  higher  collapse  load 
than  is  seen  experimentally. 

Table  5.3  shows  the  comparison  for  the  analytical 
results  from  STAGS.  Fig  5.37  shows  the  comparison  between 
12"  panels  and  8"  panels  (assuming  similar  boundary 
conditions  and  geometric  imperfections).  As  is  shown  in 
Table  5.3,  the  simply  supported  panels  were  not 
significantly  affected  by  a  reduction  in  panel  width.  It 
appears  that  the  geometric  imperfection  is  dominating  panel's 
response,  thereby  minimizing  the  effect  of  changing  the 
panel  width.  However,  the  unsupported  panels  are  affected 
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Total  Applied  Load  (lbs) 


Fig  5.31  Load-Displacement,  Simply  Supported  Panels  with 
No  Cutout 
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Fig  5.32  Load-Displacement,  Simply  Supported  Panels  with 
4"  Cutout,  Centered 
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Fig  5.33  Load-Displacement,  Simply  Supported  Panels  with 
4"  Cutout,  1"  Offset 
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Load-Displacement,  Unsupported  Panels  with 
No  Cutout 


5-49 


Total  Applied  Load 


800.00  H 


♦00.00  w 


0.00 


I  I  I  r  r  I  I  I  I  I  I  T  f  I  1  I  f  X  I  I  I  I  I  I  I  I  | 

0.01  0.01  0.02  0.02 

U  Displacement  (in) 


Fig  5.35  Load-Displacement,  Unsupported  Pane 
4"  Cutout,  Centered 
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Total  Applied  Load  (lbs) 


Fig  5.36  Load-Displacement,  Simply  Supported  Panels  with 
4"  Cutout,  1”  Offset 
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by  reducing  the  panel  width.  Combining  this  with  the 
presence  of  a  cutout  creates  a  significant  reduction  of 
total  collapse  load,  approximately  70%. 

Table  5.3  Panel  Results  by  Comparing  Panel  Geometry 

Collapse  Load,  lbs 


Panels 

12" 

8" 

Reduction 

SS,  NC 

4950 

3700 

23.5% 

SS,  CC 

2521 

2178 

13.6% 

SS,  OC 

2561 

2465 

3.7% 

US,  NC 

2260 

1350 

38.7% 

US,  CC 

1186 

323 

72.8% 

US,  OC 

1382 

460 

66.9% 

NOTE:  SS  =  Simply  Supported,  US  =  Unsupported, 

NC  =  No  Cutout,  CC  =  Centered  Cutout, 

OC  =  Off  center  Cutout 


It  should  be  noted  that  the  unsupported,  8”  panels  have 
demonstrated  the  greatest  sensitivity  in  reduction  of  total 
collapse  load  and  nonlinear  response  to  the  presence  of 
geometric  imperfections,  boundary  conditions,  panel  geometry 
or  any  combination  of  these  three  parameters.  However, 
because  of  this  high  degree  of  sensitivity,  and  combining 
these  parameters,  characterizing  the  response  of  these 
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panels  is  extremely  complex.  Further  work  needs  to  be 
pursued  before  an  adequate  analysis  can  be  accomplished 

5 . 2  Analytical  R^^sults  -  SHELL 

At  this  time,  only  one  model  has  been  successfully  run 
using  the  SHELL  program  with  the  full  model  using  a  fine 
mesh.  This  model  is  the  8"  panel  with  a 

4"  cutout,  centered  and  is  unsupported  along  the  vertical 
edges  (refer  Tables  H.3-4).  Pour  other  models,  full  scale 
but  with  coarser  meshes  (1"  elements),  were  completed. 

These  models  converged  to  solutions  from  50-80%  higher  than 
the  experimental  values  (refer  to  Table  5.4).  As  expected, 
using  a  mesh  that  modeled  half  of  the  panel,  SHELL 
converged  to  a  solution  with  a  higher  collapse  load  than 
STAGS  (from  75-90%  higher  than  STAGS).  Table  5.5  compares 
these  results.  The  solution  for  SHELL  is  quite  good 
(312.8  lbs)  when  proper  modeling  techniques  are  utilized.  T! 
load-displaceme  .t  curve  for  this  model  is  shown  in  Fig  5.38. 
At  collapse,  w  was  found  to  equal  in  magnitude  m  29.2° 

X 

or  0.51  radians.  This  would  indicate  that  through-the- 
thickness  shear  is  present  and  should  relax  the  panel's 
global  stiffness.  This  in  turn,  reduces  total  collapse  load 
and  enhance  the  bending  moment  that  occurs  in  the  panel  when 
it  nears  the  collapse  point.  The  SHELL  program  indicates 


Table  5.4  Comparing  SHELL  Models  with  Coarse  Mesh  Size 
(1"  Elements)  vs  Experimental  Results 


Collapse  Load,  lbs 


Panels 

Experimental 

SHELL 

%  Difference 

12"  US, 

NC 

1890 

2782.7 

47.2 

12"  US, 

CC 

1093 

16^5. 3 

48.7 

12"  US, 

OC 

1127 

1702.9 

51.1 

8"  US, 

NC 

1311 

2199.9 

67.8 

8"  US, 

CC 

300 

527.4 

75.8 

8"  US, 

OC 

403 

727.2 

81.2 

Table  5.5  Comparing  SHELL  Models  with  Half  Mesh  Size 
(1/2"  Elements)  vs  Experimental  Results 


Collapse  Load,  lbs 


Panels 

Experimental 

SHELL 

X  Difference 

12"  US,  NC 

1890 

2956.0 

56.4 

12"  US,  CC 

1093 

1893.1 

54.9 

8"  US,  NC 

1311 

2375.5 

81.2 

8"  US,  CC 

300 

552.9 

84.3 

NOTE!  SS=  Simple  Supports,  US  =  Unsupported,  HC  =  No 

Cutout,  CC  =  Centered  Cutout,  OC  =  Off  center  Cutout 
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on  naiiHrlw  imni 


SHELL 

Experimental 


Fig  5.38  Load-Displacement,  8"  Panel,  Unsupported  with 
4"  Cutout,  Centered 


this  with  the  full  model,  refined  mesh  solution  shown  in 
Fig  5.38  which  is  less  than  STAGS. 

The  other  models,  with  refined  meshes,  are  being 
pursued  at  this  moment,  and  will  be  presented  in  a  paper  at 
the  AIAA  conference  in  April,  1990.  From  all  indications, 
SHELL  is  correctly  predicting  the  behavior  of  the  nonlinear 
response  for  the  unsupported  panels.  Since  SHELL 
incorporates  through-the-thickness  shear  stress,  the 
solution  it  converges  to  is  realistic  and  accurate. 

However,  further  analysis  must  wait  until  comparisons  can  be 
made  concerning  all  of  the  STAGS,  SHELL  models  with  the 
experimental  data. 

5 . 3  Experimental  Results 

The  individual  experimental  panel  results,  as  are  the 
STAGS  and  SHELL  results,  are  tabulated  in  Appendix  H,  Tables 
H.1-5.  Load-displacement  curves  are  generated  for  all 
twelve  panels  and  are  shown  in  Figs  5.39-50.  We  can  see 
immediately  that  both  STAGS  and  SHELL  give  predictions  that 
are  stiff er  than  the  experimental  results,  which  is 
expected.  For  the  one  model  SHELL  has  converged  to,  SHELL 
predicts  a  more  flexible  panel  than  STAGS,  which  is 
expected.  This  can  be  attributed  to  the  fact  that  through 
the  thickness  shear  is  present  (due  to  the  large  angles  of 
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Fig  5.39  Load-Displacement,  12"  Panel,  Simply  Supported, 
No  Cutout 
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Total  Applied  Load  (lbs) 


3000.CX) 


STAGS 

Experimental 


Fig  5.40  Load-Displacement,  12"  Panel,  Simply  Supported, 
4"  Cutout,  Centered 
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Total  Applied  Load  (lbs) 


Fig  5.41  Load-Displacement,  12”  Panel,  Simply  Supported, 
4”  Cutcit,  1"  Offset 
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Pig  5.42  Load-Displacement,  8"  Panel,  Simply  Supported, 
No  Cutout 
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Total  Applied  Load  (lbs) 
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Fig  5.43  Load-Displacement,  8"  Panel,  Simply  Supported, 
4”  Cutout,  Centered 
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lotal  Applied  Load 


Total  Applied  Load  (lbs) 


Fig  5.45  Load-Displacement,  12"  Panel,  Unsupported, 
No  Cutout 
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Load-Displacement,  12"  Panel,  Unsupported, 
4"  Cutout,  Centered 
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Total  Applied  Load 


Fig  5.47  Load-Displacement,  12"  Panel,  Unsupported, 
4"  Cutout,  1"  Offset 
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Fig  5.48  Load-Displacement,  8"  Panel,  Unsupported, 
No  Cutout 
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Fig  5.49  Load-Displacement,  8"  Panel, 
4"  Cutout,  Centered 
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Fig  5.50  Load-Displacement,  8”  Panel,  Unsupported, 
4"  Cutout,  1”  Offset 
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rotation,  w  )  at  collapse,  which  increases  the  panel's 

X 

flexibility  (or  decreases  panel  stiffness). 

Table  5.6  shows  the  comparison  of  experimental  results 
with  the  analytical  results  from  STAGS.  Figs  5.51-54  show 
the  comparison  of  STAGS/experimental  collapse  ratios  for 
all  twelve  panels  considered.  It  is  evident  from  the 
load-displacement  curves  and  the  comparison  figures  that  the 
1986  version  of  STAGSC-1  with  the  corotational  method  is 
accurately  modeling  the  panels'  response  to  an 
axial -compression  load.  It  appears  that  STAGS  provides  the 
best  solutions  for  the  8"  panels  that  are  unsupported  along 
the  vertical  edges.  There  are  many  possible  reasons  for 
this  trend. 

One  possible  answer  can  be  shown  through  another  model . 
It  should  be  noted  that  STAGS  failed  to  converge  on  the  12" 
unsupported  panel  that  had  no  cutout  using  the  QUAF  411 
element.  At  first  a  mesh  of  1/2"  x  1/2"  was  tried.  Then  a 
coarser  mesh  of  1"  x  1"  was  tried.  Both  models  diverged  at 
w  1350  lbs,  far  short  of  the  experimental  results  for  this 
particular  panel.  Another  model  using  the  coarser  mesh  was 
tried,  but  incorporating  the  QUAF  410  element  instead.  When 
this  model  was  run,  STAGS  converged  to  a  solution  of  2298 
lbs  with  exactly  the  same  stiffness  as  the  411  element  had 


Table  5.6  Panel  Results  by  Comparing  Experiment  vs 
STAGS 


Panels 

Experiment 

STAGS 

12",  SS,  NC 

4851 

4950 

12"  SS,  CC 

2392 

2521 

12"  SS,  OC 

2406 

2561 

8",  SS,  NC 

3473 

3700 

8",  SS,  CC 

1918 

2178 

8",  SS,  OC 

1619 

2465 

12",  US,  NC 

1890 

2260 

12",  OS,  CC 

1093 

1186 

12",  US,  OC 

1127 

1284 

8",  US,  NC 

1311 

1350 

8",  US,  CC 

300 

323 

8",  US,  OC 

403 

460 

%  Difference 


+  11.9 


+  34.3 


+  16.4 


.9 


+  12.2 


+  12.4 


NOTE:  SS  =  Simply  Supported,  US  =  Unsupported,  NC  =  No 

Cutout,  CC  =  Centered  Cutout,  OC  =  Off  center  Cutout 


when  it  diverged  at  1350  lbs.  The  author  believes  the 
problem  may  lie  not  in  the  code  itself,  but  in  the  computer 
system  the  code  operates  on.  The  VAX  8550  operates  on  12-13 
significant  figures,  even  at  double  precision.  When  a  fine 
mesh  is  used,  or  an  element  with  many  DOF's  (such  as  the 
411)  the  global  stiffness  arrays  can  grow  to  several 
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thousand  equations.  As  these  equations  are  iterated  upon 
using  the  modified  Newton-Raphson  method,  the  error 
introduced  by  VAX  will  continue  to  grow  (due  to  the  smaller 
number  of  significant  figures)  until  the  model  either 
converges  to  an  unrealistic  solution,  or  completely 
diverges.  This  was  seen  by  Janisse  [21]  when  he  used  both 
the  VAX  and  the  Cyber  to  run  STAGS.  The  only  possible  way 
to  verify  this  hypothesis  is  to  run  the  same  model,  with 
the  411  element,  on  the  Cyber  (which  has  12  significant 
digits  at  single  precision  alone)  and  determine  if  the  model 
converges  or  not. 

Based  on  this  hypothesis,  the  author  believes  that 
STAGS  may  be  converging  to  accurate  solutions  on  the  8" 
models  for  exactly  the  same  reason.  These  models  require 
less  DOF,  therefore  the  matrices  iterated  upon  for  solution 
are  smaller,  and  the  induced  error  is  thereby  minimized. 

The  author  acknowledges  the  possibility  of  other  answers  as  to 
why  STAGS  gave  superior  answers  for  the  8"  models,  but 
believes  the  hypothesis  mentioned  above  is  viable. 
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Fig  5.51  STAOS/Experimental  Collapse  Ratios  for  12" 
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Fig  5.52  STAOS/Experimental  Collapse  Ratios  for  8" 
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Fig  5.53  STAOS/Experimental  Collapse  Ratios  for  12" 
Unsupported  Panels  Considered 
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Upon  studying  Figs  5.51-54,  it  is  clear  that  STAGS  is 
handling  large  cutouts,  eccentricity,  and  a  large  nonlinear 
response  due  to  free  vertical  edges  very  well.  The 
imperfection  model  used  on  the  solid  panels,  albeit  a  crude 
approximation  of  the  true,  initial  surface  imperfections  of 
the  panels,  allows  STAGS  to  better  predict  the  bending 
moments  that  occur  near  the  total  collapse  load  of  the 
panels.  However,  as  is  evident  in  Figs  5.39-5.50,  the 
imperfection  model  leaves  much  to  be  desired.  STAGS  is  only 
able  to  predict  a  small  portion  of  the  magnitude  of  the 
bending  moments  that  exist  when  the  panel (s)  enter  the 
nonlinear  portion  of  the  load-displacement  curve. 

In  spite  of  under  estimating  the  bending  moments,  STAGS 
predicts  quite  well  the  radial  displacement  pattern  that  is 
seen  experimentally  through  the  LVDT  measurements.  For 
brevity,  none  of  these  curves  are  included.  Unlike  the 
results  Tisler  [35]  reported,  STAGS  did  not  demonstrate  a 
softening  of  the  panel.  It  is  believed  that  this  is  due  to 
the  advancements  in  the  code  made  through  the  corotational 
(updated  Lagrangian  scheme)  method. 

The  experimental  radial  deformation  pattern,  while  very 
similar  to  that  seen  for  a  perfect  panel,  is  much  different 
that  what  is  predicted  using  the  initial  imperfection  model . 
Combining  this  with  the  fact  that  all  the  panels  had 
different,  random  initial  imperfection  patterns  and  yet 
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still  failed  in  the  same  manner,  we  conclude  that  large 
cutouts  dominate  the  panel (s)'  response  to  axial 
compression.  The  initial  imperfections  play  a  role  in  the 
response  only  in  those  panels  that  do  not  have  geometric 
imperfections . 

At  the  point  of  collapse,  for  the  simply  supported 
panels,  the  panel  is  seen  to  snap  into  a  new  form.  Although 
this  snap-through  is  more  dramatic  for  the  solid  panels, 
those  panels  with  4"  cutouts  did  exhibit  a  secondary  mode. 
This  was  seen  for  both  12“  and  8“  panels  (Figs  5.55-56).  As 
is  mentioned  in  [35],  apparently  enough  energy  is  retained 
by  the  vertical  supports  to  force  the  panel (s)  to  collapse, 
via  a  bifurcation-like  approach,  into  a  secondary  mode.  A 
possible  contribution  to  this  may  be  the  surface 
imperfection  along  the  cutout  edges  causing  energy  to 
buildup  in  a  different  manner  than  STAGS  predicts. 

The  solid  panel  with  free  vertical  edges  also  acted  as 
having  a  secondary  mode  (Figs  5.57-58).  It  is  thought  that 
the  panel,  since  there  are  no  vertical  supports  to  add 
stiffness,  essentially  acts  as  plate  and  therefore  exhibits 
the  secondary  mode  once  collapse  is  achieved. 

For  those  panels  with  free  vertical  edges  and  geometric 
imperfections,  the  panels  did  not  exhibit  any  snap-through 
phenomenon  (Figs  5.59-62).  Instead,  once  the  total  collapse 
load  is  achieved,  the  panel  continues  its  displacement 
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Fig  5.57  Picture  of  12"  Unsupported  Panel,  No  Cutout, 
Prior  to  Collapse  Snap-Through 
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Pig  5.61 


Picture  of  8"  Unsupported  Panel  with  4"  Cutout 
Prior  to  Collapse 
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Pig  5.62  Pictures  of  8"  Unsupported  Panel  with  4"  Cutout 
After  to  Collapse 
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bending  moments  through  the  nonlinear  portions  of  the 
pattern  until  the  loading  is  removed. 

In  all,  STAGS  and  SHELL  predicts  the  general 
displacement  pattern,  both  axially  and  radially,  very  well. 
When  the  imperfection  models  are  included,  STAGS  does 
predict  an  increase  in  load-displacement  curves.  STAGS  and 
SHELL  also  accurately  predicts,  the  effects  of  boundary 
conditions,  geometric  imperfections,  eccentricity,  and  panel 
geometry  upon  composite  shell  response  to  axial  compression. 

5 . 4  Residual  Stress 

The  results  of  the  residual  stress  test  showed  little 
change  in  the  strain  in  the  axial  direction  (w  lO  /i  strain). 
The  circumferential  direction,  however,  saw  a  relatively 
large  change  in  surface  strain  (w  75  n  strain).  How  much  of 
this  could  be  attributed  to  experimental  error  is  unclear. 
Although  the  panel  was  measured  before  cutting  while  freely 
resting,  the  curved  panel  is  very  sensitive  to  bending 
stresses.  Some  bending  stress  may  have  been  caused  by 
merely  resting  the  panel  on  the  wooden  cutout  form  (refer 
Section  3.3,  Fig  3.18). 

Assuming  the  experimental  findings  are  correct,  it 
remains  to  compare  these  stresses  to  stresses  during  axial 
compression.  The  change  in  strain  appears  to  be  in  bending 


stress  present  due  to  curing  the  panel  during  the 
manufacturing  process.  This  means  that  the  maximum  change 
in  stress  occurred  at  the  outer  surface  where  the  strain 
rosettes  were  positioned.  This  change  in  bending  stress  was 
of  equal  and  opposite  magnitude  on  the  opposing  side  of  the 
panel,  and  zero  at  the  center  of  the  panel. 

Since  the  change  in  stress  in  the  axial  direction  was 
negligible  («  14  psi),  interest  is  in  the  stress  in  the 
circumferential  direction.  The  change  in  stress  in  the 
circ\imf erential  direction  (at  the  outermost  ply)  that 
corresponds  to  the  measured  change  in  strain  varied  from  450 
psi  at  the  corners  to  750  psi  at  the  center  of  the  cutout. 

A  plot  of  the  circumferential  stress  at  collapse  is  shown  in 
Fig  5.63,  where  the  stress  magnitude  is  measured  for  the 
cutout  edge,  positive  outward.  It  is  seen  that  the  change 
in  stresses  are  much  smaller  in  magnitude  than  the  stresses 
from  loading.  These  results  differ  greatly  from  Tisler's 
[31]  results.  The  author  believes  this  is  due  mainly  to  the 
experimental  approach  taken  between  this  effort  and 
Tisler's.  Tisler  used  a  single  strain  rosette,  placed  in 
the  center  of  the  cutout.  This  single  rosette  could,  at 
best,  give  a  rough  global  approximation  of  the  change  in 
stress  due  to  removing  the  cutout.  The  approach  used  in 
this  effort  incorporated  16  strain  rosettes  (Fig  5.64)  on 
the  front  and  back  of  the  panel.  These  rosettes  were  placed 
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0.5"  away  from  the  cutout  to  allow  clearance  for  the  router 
used  in  removing  the  cutout  (refer  Section  3.3,  Fig  3.19). 
Although  this  is  relatively  small  distance,  it  may  be  great 
enough  to  begin  relieving  the  change  in  stresses.  However, 
the  author  believes  the  results  generated  from  this 
experimental  approach  give  a  clearer  picture  as  to  the 
magnitude  and  shape  of  the  residual  stress  than  before.  It 
is  seen  to  be  at  least  50%  less  than  stresses  from  loading 
at  collapse. 

At  this  point  it  is  unclear  as  to  the  importance  this 
change  in  stress  is  to  the  panel(s)'  collapse.  How  much  (if 
any)  of  the  error  in  the  stags  and  SHELL  solutions  can  be 
attributed  to  residual  stress  changes  is  unclear.  We  have 
seen  that  based  on  the  previous  work,  along  with  this 
effort,  this  is  a  subject  that  warrants  further  study.  It 
is  evident  that  residual  stress  changes  exist.  What  needs 
to  be  determined  is  the  role  these  stress  changes  play  in 
the  collapse  of  cylindrical  composite  shells.  A  statistical 
study  should  be  made  in  determining  the  stress 
concentrations  along  the  edges  of  the  cutout  (as  was  done  in 
this  effort)  should  provide  a  clearer  picture  as  to  the 
importance  of  considering  residual  stress  in  the  collapse 
analysis.  The  the  only  references  available  on  this  subject 
are  for  isotropic  materials  with  circular  holes  in  plates 
and  do  not  consider  cylindrical  composite  shells  with  large 
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cutouts.  Since  this  is  the  case,  a  methodology  of 
incorporating  these  residual  stresses  into  the  collapse 
analysis  along  with  a  larger  data  base  need  to  be  developed. 

One  final  note,  there  is  some  concern  about  inducing 
residual  stress  by  forcing  the  experimental  panels  to  fit 
the  curvature  of  the  test  fixture.  While  the  panel(s)  are 
curing  in  the  autoclave,  the  original  radius  of  curvature 
is  no  longer  held  constant  due  to  the  anisotropy  of  the 
laminate.  Hence  the  panel,  once  curing  was  completed,  now 
has  a  random  radius  of  curvature  throughout.  However,  the 
test  fixture  was  machined  with  a  constant  radius  of 
curvature  of  12"  (the  original  curvature  for  the  laminate 
when  it  was  mar.^«-\*ctured) .  By  forcing  the  panel  into  this 
constant  radias  of  curvature,  it  is  felt  that  internal 
stresses  may  be  induced,  thereby  decreasing  panel  stiffness 
and  ac  the  same  time  increasing  the  bending  moment  seen  in 
the  nonlinear  panel  response.  By  referring  to  Section  3.3, 
Fig  3.17,  the  reader  is  reminded  that  a  reading  was  taken 
from  the  strain  rosettes  while  the  panel  was  strapped  to  the 
hole  cutting  support.  It  is  felt  that  this  closely 
parallels  the  situation  when  the  panel  was  fit  into  the 
experimental  test  fixture.  An  analysis  was  accomplished, 
similar  to  the  actual  residual  stress  test  and  is  shown  in 
Fig  5.65.  It  was  determined  that  forcing  the  panel  to 
assume  a  constant  radius  of  curvature  generated  a  maximum  of 
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250  lbs.  Granted,  this  is  a  small  amount  (only  9%  of  the 
total  collapse  load)  but  it  is  of  enough  magnitude  that  it 
can  be  considered  in  the  analysis  if  so  desired. 


6 .  CONCLUSIONS 


1.  It  is  still  unclear  whether  the  vertical  boundary 
condition  of  circumferential  movement  (v)  is  fixed  or 
free.  From  the  analysis  presented  in  Chapter  5,  it 
appears  the  true  boundary  condition  for  the 
experimental  panels  is  somewhere  in  between  the  two 
extremes.  The  v  free  condition  does  not  reduce  panel 
stiffness,  but  it  will  reduce  the  total  collapse  load 
for  the  panels.  For  those  panels  with  no  cutouts  the 
reduction  is  w  8%.  For  those  panels  with  4“  cutouts, 
the  reduction  is  18%.  This  is  because  the  majority  of 
the  load  is  carried  along  the  vertical  edge  if  there  is 
a  large  cutout  present.  With  a  solid  panel,  the  load 
remains  evenly  distributed. 

2.  Small  radial  imperfections  effects  become  less  crucial 
when  cutout  are  present,  and  are  almost  negligible  when 
large  cutouts  are  introduced.  The  impact  of  these 
imperfections  is  minimal  when  compared  to  the  effect  of 
the  cutouts  on  the  panel(s)*  response.  However,  these 
small  imperfections  do  play  a  role  in  those  panels  no 
having  a  cutout.  STAGS  has  shown  that  using  an  initial 
imperfection  model  using  three  half  sine  waves  in  both 
the  radial  and  circumferential  directions  with  an 
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amplitude  of  one  ply  thickness  (0.005“  )  indicates  a 
reduction  in  stiffness  for  both  simply  supported  panels 
(5-10%)  and  unsupported  panels  (15-35%).  For  the 
unsupported  panels,  this  reduction  in  stiffness 
included  a  marked  increase  in  bending  for  the  nonlinear 
portion  of  the  load-displacement  curves. 

The  presence  of  geometric  imperfections  (cutouts) 
reduces  the  panel’s  total  collapse  load  by  m  50%  (up  to 
77%  for  the  unsupported  8“  panels).  Offsetting  the 
cutout  to  introduce  eccentricity  did  not  reduce  the 
collapse  load  as  expected.  Instead,  STAGS  predicted  a 
slight  increase  in  total  collapse  load  (3%  for  simply 
supported  panels  and  7-10%  for  unsupported  panels) 
which  is  verified  by  experimental  results.  This 
increase  in  loading  may  be  due  to  the  redistribution  of 
stress  that  occurs  between  the  cutout  and  the  vertical 
edges.  If  the  reader  considers  the  12“  simply 
supported  panel  with  a  4“  centered  cutout  as  the 
baseline:  then  as  the  boundary  is  moved  closer,  the 
stress  (N  )  redistribution  occurs  at  a  lower  percentage 

X 

of  the  total  collapse  load  than  is  seen  for  the 
baseline.  Thus  in  the  case  of  eccentricity,  the 
smaller  “column“  or  side  of  the  panel  redistributes  its 
stress  at  a  lower  magnitude  of  stress,  and  the  wider 


column  will  keep  an  even  distribution  at  a  higher 
magnitude  of  stress.  There  may  be  some  localized 
buckling  due  to  the  smaller  column  reaching  its  maximum 
stress  at  a  lower  magnitude  of  axial  loading,  but  from 
a  global  viewpoint,  the  panel  will  respond  with  an 
essentially  increased  stiffness. 

4.  The  vertical  boundary  conditions  have  a  significant 
impact  on  the  total  collapse  load.  By  removing  the 
vertical  supports,  total  collapse  load  of  the  panel 
(ignoring  the  effect  of  cutouts  or  panel  width)  is 
reduced  by  at  least  40%  and  up  to  75%.  This  is  due  to 
essentially  removing  the  ±  45°  plies  from  the  panel’s 
stiffness  matrix  (which  is  half  of  the  ply  lay-up) 
since  there  is  no  vertical  support  to  constrain  these 
plies.  Hence,  the  *»  50%  reduction  in  total  collapse 
load  for  all  six  panels  when  the  vertical  supportes  are 
removed . 

5.  By  reducing  the  panel’s  width  and  thereby  placing  the 
effects  of  the  vertical  boundary  condition  near  the 
cutout,  a  reduction  in  collapse  load  is  seen.  The 
effect  of  panel  geometry  is  not  as  significant  to  the 
panel’s  response  when  considering  cutouts  or  boundary 
conditions  of  the  panel  for  the  simply  supported 
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panels.  This  is  not  true  for  the  unsupported  panels. 

For  the  simply  supported  panels,  geometric 
imperfections  (cutouts)  and  removing  the  boundary 
conditions  have  a  greater  impact  than  reducing  the 
Midth  of  the  panel .  The  unsupported  panels  are  very 
sensitive  to  the  presence  of  cutouts  and  reducing 
panel  width.  Both  conditions  have  approi.^mately  the 
same  result,  a  large  decrease  in  panel  stiffness  and 
therefore  a  significant  decrease  in  total  collapse 
load . 

The  unsupported  panels  demonstrated  the  greatest 
sensitivity  to  presence  of  cutouts.  By  removing 
one-third  to  half  of  the  effective  cross-section  of  the 
panel  to  carry  a  load,  the  panel  could  see  a  reduction 
in  total  collapse  load  of  up  to  7S\.  The  same  holds 
true  from  reducing  panel  width  from  12“  to  8“.  The 
response  of  these  panels  is  highly  nonlinear  with 
significant  rotations  occurring  at  the  point  of  collapse 
along  the  free  vertical  edge.  These  large  rotations 
may  indicate  the  presence  of  through-the-  thickness 
shear  which  reduces  panel  stiffness  and  therefore  total 


collapse  load. 
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STAGS  appears  quite  able  to  handle  large  displacements 
and  rotations  exhibited  by  the  panels.  Although  STAGS 
may  predict  higher  stiffnesses  for  the  panels,  this  is 
to  be  expected.  The  imperfection  model  used  forces  the 
random  imperfections  known  to  exist  on  the  experimental 
panels  into  a  uniform  pattern  through  the  sine  wave 
approximation  which  minimizes  the  amount  of  bending 
moment  introduced  into  the  model  as  it  nears  the 
collapse  load.  Residual  stress  (due  to  curing  and 
enforcing  a  constant  radius  of  curvature)  is  not 
considered.  For  the  simply  supported  panels,  the  v 
free  or  fixed  boundary  condition  is  still  unclear.  Yet 
STAGS  accurately  predicts  the  panel’s  response  in  terms 
of  load,  axial  displacement,  and  radial  displacement. 

8.  SHELL  predicts  more  accurately  than  STAGS  for  those 
panels  that  undergo  large  displacements  (m  .02")  and/or 
large  rotations  (w  15-17®).  This  is  due  to  transverse 
shear  reducing  the  panel’s  stiffness.  This  reduction 
in  stiffness  is  seen  through  a  reduction  in  total 
collapse  load,  an  increase  in  the  bending  moments,  and 
an  increase  in  the  magnitude  of  radial  displacments . 

9.  Both  STAGS  and  SHELL  are  useful  analysis  tools.  STAGS 
is  more  advantageous  for  analyzing  general  shell  shapes 
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10. 


(either  isotropic  or  anisotropic  material)  that  undergo 
small,  moderate,  and  even  large  displacements  and/or 
rotations.  On  the  other  hand,  STAGS  ignores 
three-dimensional  effects  and  cannot  adequately  conduct 
nonlinear  vibrational  analysis.  SHELL  has  the 
advantage  of  incorporating  transverse  shear  effects 
into  its  analysis,  its  strain-displacement 
relationships  are  based  upon  the  exact  nonlinear,  shell 
strain-displacement  relationships,  and  the  additional 
capability  of  nonlinear  vibrational  analysis  (both 
static  and  dynamic).  However,  SHELL  can  only  analyze 
cylindrical  shell  structures  or  plate  structures. 

Residual  stress  is  created  by  removing  the  cutout.  It 
is  fairly  uniform  in  the  circumferential  direction  and 
practically  nonexistent  in  the  axial  direction.  The 
magnitude  of  this  stress  is  approximately  750  lbs 
(assuming  worst  case).  When  comparing  this  to  the 
total  collapse  load  (2521  lbs)  it  is  evident  that 
residual  stress  does  not  dominate  nor  greatly  influence 
the  panel’s  response.  It  is  felt  the  presence  of  the 
cutout  itself  dominates  the  reduction  in  total  collapse 
load  that  is  seen.  However,  residual  stress  should  not 
be  ignored. 
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7 .  RECOMMENDATIONS 


1.  Further  analysis  need  to  be  accomplished  on  the  initial 
imperfection  model.  It  is  believed  the  additional 
bending  seen  in  the  nonlinear  resonse  of  the  panel  is 
due  to  radial  surface  imperfections.  One  method  is 
mark  out  the  finite  element  model  on  the  surface  of  the 
panel .  Assume  a  uniform  radius  of  curvature  through 
the  midplane  of  the  panel.  Then  determine,  at  each 
node  point  on  the  panel  the  change,  in  terms  of 
displacement,  from  this  uniform  plane  of  constant 
radius  of  curvature.  Then  input  these  displacements 
into  the  model  used  in  STAGS  as  initial  imperfection 
and  the  solve  the  model  for  equlibrium  as  load  is 
increase.  This  should  provide  a  solution  that 
approximates  the  experimental  results  more  closely. 
Another  method  is  to  incorporate  the  techniques  Starnes 
[30]  used  into  STAGS.  The  author  realizes  this  is  not 
a  simple  problem,  but  it  is  felt  to  be  important  enough 
to  be  pursued. 

2.  More  analysis  needs  to  be  accomplished  on  the 
cylindrical  shells  that  have  free  vertical  edges. 

SHELL  indicates  that  through-the-thickness  shear  is 
important  at  the  onset  of  collapse.  This  shear  effect 
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could  reduced  panel  stiffness  and  thereby  decrease 
total  collapse  load.  It  may  also  enhance  the  bending 
moment  that  begins  to  dominate  the  panel’s  response  in 
the  nonlinear  region. 

Verification  is  still  needed  on  the  experimental 
vertical  bounary  conditions,  whether  v  is  free  or 
fixed.  One  solution  is  model  the  panel  such  that  on 
the  vertical  edges,  every  other  node  has  v  free  and  the 
rest  have  v  fixed.  The  author  acknowledges  this  is  not 
the  actual  circumstances  seen  in  the  experimental 
setup,  but  this  model  could  give  a  better  approximation 
for  the  true  boundary  conditions. 


8 .  BIBLIOGRAPHY 


1-  Ugral ,  A.  C.,  Stresses  in  Plates  and  Shells . 

McGraw-Hill  Book  Co.,  1981. 

2 .  Saada ,  A .  S . ,  Elasticity  Theory  and 
Applications .  Pergommen  Press,  1974. 

3.  Sanders,  J.  L.  "Nonlinear  Theories  for  Thin  Shells," 

Qtr .  of  Applied  Math.  Vol  XXI,  No.  1  (1962). 

4.  Leissa,  A.  U.  Buckling  of  Laminated  Composite  Plates 

and  Shell  Panels .  AFWAL-TR-85-3069 .  Air  Force  Flight 

Dynamics  Laboratory,  Uright-Patterson  AFB  OH,  June 
1985. 

5.  Sobel ,  L.  H.,  Weller,  T.,  and  Agarwal ,  B.  L.,  "Buckling 
of  Cylindrical  Panels  Under  Axial  Compression," 

Computers  and  Structures.  6:  29-35  (February  1976). 

6.  Almroth,  B.  0.  "Influence  of  Edge  Conditions  on  the 
Stability  of  Axially  Compressed  Cylindrical  Shells," 

AIAA  Journal .  4:  134-14  (1969). 

7.  Andreev,  L.  V.  et  al . ,  "Nonlinear  Deformation  of  a 
Cylindrical  Panel  Under  Axial  Compression,  Soviet 
Applied  Mechanics.  17:  270-274  (September  1981). 

8.  Rehfield,  L.  U.  and  Hallauer,  U.  L.,  Jr.,  "Edge 
Restraint  Effect  on  Buckling  of  Compressed  Curved 
Panels,"  AIAA  Journal .  6:  187-189  (January  1980). 

9.  Becker,  M.  L.,  Palazotto,  A.  N.,  and  Khot,  N.  S., 
"Instability  of  Composite  Panels,"  Journal  of  Aircraft, 
18:  739-743  (September  1981). 

10.  Tvergaard,  V.,  "Buckling  of  Elastic-Plastic  Cylindrical 
Panel  Under  Axial  Compression,"  International  Journal 
of  Solid  Structures .  13:  957-970  (1977). 

11.  Brogan,  F.  and  Almroth,  B.  0.,  "Buckling  of  Cylinders 
with  Cutouts,"  AIAA  Journal .  8:  236-238  (February  1970). 

12.  Garashchuk,  I.  N.,  and  Chernyshenko,  I.  S.,  "Numerical 
Investigation  of  the  Stress  State  of  a  Cylindrical 
Shell  with  a  Circular  Hole,"  Soviet  Applied  Mechanics . 
15:  484-488  (December  1979). 


13.  Almroth,  B.  0.  and  Holmes,  A.  M.  C.,  “Buckling  of 
Shells  with  Cutouts  Experiment  and  Analysis,” 
International  Journal  of  Solid  Structures .  8:  1057-1071 
( August  1972 ) . 

14.  Almroth,  B.  0.,  et  al . ,  “Stability  of  Cylinders  with 
Circular  Cutouts,”  AIAA  Journal .  11:  1582-1584 

( November  1973  ) . 

15.  Almorth,  B.  0.,  and  Brogan,  F.  A.  “Bifurcation  Buckling 
as  an  Approximation  of  the  Collapse  Load  for  General 
Shells,”  AIAA  Journal .  10:  463-467  (April  1972). 

16.  Becker,  M.  L.,  Palazotto,  A.  N.,  and  Khot ,  N.  S., 
“Experimental  Investigation  of  the  Instability  of 
Composite  Cylindrical  Panels,”  Experimental  Mechanics . 
22:  372-376  (October  1982). 

17.  Bauld,  N.  R.,  Jr.,  and  Satyamurthy,  K.,  Collapse  Load 

Analysis  for  Plates  and  Shells .  AFFDL-TR-79-3038 .  Air 

Force  Flight  Dynamics  Laboratory,  Wright-Patterson  AFB 
OH ,  December  1970 . 

18.  Harper,  James  G.  Buckling  Analysis  of  Laminated 
Composite  Circular  Cylindrical  Shells .  MS  Thesis, 
AFIT/GAE/AA/78D-8 .  School  of  Engineering,  Air  Force 
Institute  of  Technology  (AU),  Wright-Patterson  AFB  OH, 
December  1979 . 

19.  Becker,  Marvin  L.  Analytical/Experimental 
Investigation  of  the  Instability  of  Composite 
Cylindrical  Panels.  MS  Thesis,  AFIT/GAE/AA/79D-3 . 
School  of  Engineering,  Air  Force  Institute  of 
Technology  (AU),  Wright-Patterson  AFB  OH,  December 
1979. 

20.  Hebert,  John  S.  Analytical/Experimental  Linear 
Bifurcation  of  Curved  Cylindrical  Composite  Panels .  MS 
Thesis,  AFIT/GAE/AA/82D-14 .  School  of  Engineering,  Air 
Force  Institute  of  Technology  (AU),  Wright-Patterson 
AFB  OH,  December  1982. 

21.  Janisse,  T.  C.  and  Palazotto,  A.  N..  “Collapse  Analysis 
of  Cylindrical  Composite  Panel  with  Cutouts,”  AIAA 
Journal  of  Aircraft,  21:  731-733  (September  1984). 

22.  Lee,  C.  E.,  and  Palazotto,  A.  N.,  “Nonlinear  Collapse 
Analysis  of  Composite  Cylindrical  Panels  with  Small 
Cutouts  or  Notches,"  Composite  Structures .  4:  217-229 
( 1985 ) . 


8-2 


23.  Hermsen,  M.  F.,  and  Palazotto,  A.  N.,  "The  Effects  of 
Cutout  Location  and  Material  Degradation  on  the 
Collapse  of  Composite  Cylindrical  Panels,"  Non-Linear 
Analysis  and  NDE  of  Composite  Material  Vessels  and 
Components .  (Edited  by  Hui ,  D.,  Duke,  J.  C.,  and  Chung, 
H.)  ASME  PVP  vol.  115  and  NDE  vol .  3:  43-57  (1986). 

24.  Dennis,  S.  T.  and  Palazotto,  A.  N.,  "The  Effects  of 
Large  Movement  on  a  Composite  Cylindrical  Shell  with 
Cutouts,"  30th  AIAA  Conf . .  Mobile  AL ,  paper  no. 89-1398, 
to  appear  in  Nonlinear  Mechanics  (1989). 


25.  Linneman,  P.  E.  Vibration  and  Buckling  Characteristics 
•  of  Composite  Cylindrical  Panels  Incorporating  the 

Effects  of  a  Higher  Order  Shear  Theory .  MS  Thesis, 
AFIT/GA/AA/88D-6 .  School  of  Engineering,  Air  Force 
Institute  of  Technology  (AU),  Wr ight-Patterson  AFB  OH, 
December  1988 . 


•  26.  Cook,  R.  D.  Concepts  and  Applications  of  Finite 

Element  Analysis  for  Plates  and  Panels .  (Second 
Edition).  John  Uiley  and  Sons,  1981. 


27. 


User  *s  Manual  for  ST  AGS .  Volume  1  Theory.  Structural 
Mechanics  Laboratory,  Lockheed  Palo  Alto  Research 
Laboratory,  Palo  Alto  CA,  1981. 


28.  Zienkiewicz,  0.  C.  The  Finite  Element  Method  (The 
third,  expanded  and  revised  edition  of  The  Finite 
Element  Method  in  Engineering  Sciences).  McGraw-Hill 
Book  Company  (UK)  Limited,  1977. 

29.  Lee,  Cathy  E.  Numerical  Determination  of  the  Effects 
of  Boundary  Conditions  on  the  Instability  of  Composite 
Panels  with  Cutouts .  MS  Thesis,  AFIT/GAE/AA/83D-4 . 
School  of  Engineering,  Air  Force  Institute  of 
Technology  (AU),  Ur ight-Patterson  AFB  OH,  December 
1983. 

30.  Starnes,  J.  H.  and  Knight,  N.  F.,  "Postbuckling 
Behavior  of  Selected  Graphite-Epoxy  Cylindrical  Panels 
Loaded  in  Axial  Compression,”  27th  Structures, 
Structural  Dynamics,  and  Materials  Conference,  San 
Antonio  TX ,  May  19-21  1986,  AIAA  86-0881. 

31.  Knight,  N.  F.,  and  Starnes,  J.  H.,  "Postbuckling 
Behavior  of  Axially  Compressed  Graphite  Epoxy 
Cylindrical  Panels  with  Circular  Holes,"  Collapse 
Analysis  of  Structures .  (ed.  by  L.  H.  Sobel )  ASME,  PVP 
Vol .  84 ,  ( 1984  )  . 


8-3 


•  32.  Palazotto,  A.  N.,  and  Tisler,  T.  W.,  “Experimental 

Collapse  Determination  of  Cylindrical  Composite  Panels 
with  l^Large  Cutouts  under  Axial  Load,"  Composite 
Structures .  4:  61-78  (1989). 

33.  Palazotto  A.  N.,  and  Tisler.  T.  UJ.,  "Considerations  of 

•  Cutouts  in  Composite  Cylindrical  Panels,"  Computers  & 

Structures .  Vol .  29,  No.  6,  pp  1101-1110,  (1988). 

34.  Tisler,  T.  U.,  and  Palazotto,  A.  N.,  "Insight  into  the 
Collapse  of  Composite  Cylindrical  Panels  with  Cutouts: 
Analysis  and  Experimentation,"  Recent  Advances  in 

•  Structural  Dynamics .  (ed.  by  H.  Chung  and  H.  D.  Fisher) 

ASME ,  PVP  Vol .  124  ( 1987  )  . 

35.  Tisler,  Thomas  W.  Collapse  Analysis  of  Cylindrical 
Composite  Panels  with  Large  Cutouts  Under  an  Axial 
Load.  MS  Thesis,  AFIT/GAE/AA/86D-1 .  School  of 

•  Engineering,  Air  Force  Institute  of  Technology  (AU), 

Wright-Patterson  AFB  OH,  December  1986. 

36.  Egan,  G.  S.,  and  Palazotto  A.  N.,  "The  Analysis  of  a 
Composite  Shell  Structure,"  30th  AIAA/ASME  Structures, 
Structural  Dynamics  and  Materials  Conference,  Mobile 

•  AL,  Paper  No.  89-1297,  to  appear  in  AIAA  Journal 

( 1989). 

37 .  Egan ,  Gregory  S .  Nonlinear  Finite  Element  Analysis  of 
a  General  Composite  Shell .  MS  Thesis, 
AFIT/GAE/AA/88D-12 .  School  of  Engineering,  Air  Force 

•  Institute  of  Technology  (AU),  Wright-Patterson  AFB  01^ 

December  1988. 

38.  Palazotto,  A.  N.,  "An  Experimental  Study  of  a  Curved 
Composite  Cutout , "  Composite  Materials:  Testing  and 
Design  (Eigth  Conference),  ASTM  STP  972.  Ed.  by  J.  D. 

•  Whitcomb,  pp  191-202  (1988). 

39.  Riks,  E.,  "On  the  Numerical  Solution  of  Snapping 
Problems  in  the  Theory  of  Elastic  Stability,"  SUDAAR 
No.  401,  Stanford  University,  Stanford,  CA,  1970. 

•  40.  Dennis,  S.  T.,  and  Palazotto  A.  N.,  "The  Effects  of 

Large  Movement  on  a  Composite  Cylindrical  Shell  with 
Cutouts",  30th  AIAA/ASME  Structures,  Structural 
Dynamics  and  Materials  Conference,  Mobile  AL,  Paper  No. 
89-1398,  AIAA  Journal  (Oct  1989). 


8-4 


41.  Dennis,  Scott  T.  Large  Displacement  and  Rotational 
•  Formulation  for  Laminated  Cylindrical  Shells  Including 

Parabolic  Transverse  Shear .  Ph.D.  Dissertation, 
AFIT/DS/AA/88-1 .  School  of  Engineering,  Air  Force 
Institute  of  Technology  (AU),  Wright-Patterson  AFB  OH, 
May,  1988  ( AD-A194871 ) . 


•  42.  Sobel ,  L.  H.  and  Thomas,  K.  Evaluation  of  the  STAGSC-1 

Shell  Analysis  Program .  WARD- 1 088 1 .  Westinghouse 
Advanced  Reactors  Division,  Madison  PA,  August  1981. 


43.  Almroth,  B.  0.  and  Brogran,  F.  A.  Numerical  Procedures 
for  Analysis  of  Structural  Shells .  AFWAL-TR-80-3128 . 

#  Air  Force  Flight  Dynamics  Laboratory,  Wright-Patterson 

AFB  OH,  March  1981. 


44.  Brush,  D.  0.,  and  Almroth,  B.  O.  Buckling  of  Bars . 
Plates .  and  Shells .  McGraw-Hill,  1975. 

•  45.  Jones,  R.  M.  Mechanics  of  Composite  Materials. 

Hemisphere  Publishing  Corporation,  1975. 


46.  Rankin,  C.  C.,  ^  al-  "Enhancements  to  the  STAGS 
Computer  Code."  NASA  Contracotr  Report  4000,  November 
1986. 

47.  Rankin,  C.  C.  and  Brogan,  F.  A.  "An 
Element-Independent  Corotational  Procedure  for  the 
Treatment  of  Large  Rotations,"  in  Collapse  Analysis  of 
Structures .  ed.  by  L.  H.  Sobel,  ASME,  PVP  Vol .  84, 

1984. 

48.  Bathe,  K.  J.  and  Ho,  L.  W.  "Some  Results  in  the 
Analysis  of  Thin  Shell  Structures,"  Nonlinear  Finite 
Element  Analysis  in  Structural  Mechanics .  Ed.  by  W. 
Wunderlich,  E.  Stein,  and  K.  J.  Bathe,  1981,  Springer  - 
Verlag  Berlin,  Heidelberg,  NY,  (125-150). 

49.  Riks,  E.  "An  Incremental  Approach  to  the  Solution  of 

Snapping  and  Buckling  Problems,"  International  Journal 
of  Solids  and  Structures .  15:  529-551  (1979). 


50.  Wempner ,  G.  A.  "Discrete  Approximations  Related  to 

Nonlinear  Theories  of  Solids,"  International  Journal  of 
Solids  and  Structures.  7'  1581-1599  (1971). 


8-5 


51.  Tsai,  C.  T.  and  Palazotto,  A.  N.  “A  Modified  Riks 
Approach  to  Composite  Shell  Snapping  Using  High-Order 
Shear  Deformation  Theory,"  accepted  for  publication  in 
Journal  of  Computers  and  Structures  ( 1989 ) . 

52.  Horban,  B.  A.  The  Effects  of  Through-the-Thickness 
Delaminations  on  Curved  Composite  Panels.  MS  Thesis, 
AFIT/GAE/AA/85D-8 .  School  of  Engineering,  Air  Force 
Institute  of  Technology  (AU),  Wright-Patterson  AFB  OH, 
Dec  1985 . 

53.  Wilder,  B.  L.  A  Study  of  Damage  Tolerance  in  Curved 
Composite  Panels.  MS  Thesis,  AFIT/GAE/AA/88D-2 . 

School  of  Engineering,  Air  Force  Institute  of 
Technology  ( AU  )  ,  Wright-Patterson  AFB  OH,  Dec  1988. 

54.  Broek,  D.  Elementary  Engineering  Fracture  Mechanics. 
Martinus  Nijhoff  Publishers,  Dordrecht,  Netherlands, 
1987. 

55.  Crisfield,  M.  A.,  “A  Fast  Incremental/Iterative 
Solution  Procedure  That  Handles  ’Snap-Through’," 
Computers  and  Structures .  13:  55-62  (1981). 

56.  Bowlus,  J.  A.,  Palazotto,  A.  N.,  and  Whitney,  J.  M., 
"Vibration  of  Symmetrically  Laminated  Rectangular 
Plates  Considering  Deformation  and  Rotary  Inertia," 

AIAA  Journal .  25:  1500-1511  (Nov  1987). 

57.  Bowlus,  John  A.  The  Determination  of  the  Natural 
Frequencies  and  Mode  Shapes  For  Anisotropic  Laminated 
Plates  Including  the  Effects  of  Shear  Deformation  and 
Rotary  Inertia .  MS  Thesis,  AFIT/GA/AA/85S-1 .  School 
of  Engineering,  Air  Force  Institute  of  Technology  (AU), 
Wright-Patterson  AFB  OH,  September  1985. 

58.  Palardy,  Real  F.  The  Buckling  and  Vibration  of 
Composite  Plates  Using  the  Levy  Method  Considering 
Shear  Deformation  and  Rotary  Inertia .  MS  Thesis, 
AFIT/GAE/AA/87D-16 .  School  of  Engineering,  Air  Force 
Institute  of  Technology  (AU),  Wright-Patterson  AFB,  OH, 
December  1987 . 

59.  Reddy,  J.  N.,  and  Liu,  C.  F.,  "A  Higher-Order  Shear  and 
Deformation  Theory  of  Laminated  Elastic  Shells," 

Int  J  Ena  Sci.  23:  319-330  (Mar  1985). 


8-6 


Appendix  A:  Derivation  of  STAGSC-1  Nonlinear 
Strain  Displacement  Equations 


Since  one  of  the  primary  strengths  of  the  STAGS 
computer  programs  is  its  ability  to  analyze  shells 
nonlinearly,  a  derivation  of  the  nonlinear  kinematics  is 
presented . 

To  start  the  derivation,  a  few  assumptions  need  to  be 
made.  STAGS  uses  plate  elements  to  model  a  shell  structure. 
These  plate  elements  are  considered  thin  so  that  the 
in-plane  displacements,  u  and  v,  and  the  normal 
displacements,  w,  are  functions  of  only  two  space  variables 
[27] .  A  plane  stress  assumption  is  made  whereby  Y  »  y  , 

xy  yz 

e  ,  and  o  are  considered  neglible.  And  finally,  the 
2  2 

Kirchhoff-Love  hypothesis  applies  for  strains  away  from  the 
mid-plane . 

Consider  a  line  segment  as  shown  in  Fig  A.l  oriented  in 
its  undeformed  and  deformed  state  (represented  by  an  4:). 

From  Fig  A.l,  it  is  evident  that  the  deformed  state  can  be 
related  to  the  undeformed  state  by  [27] : 

X  =  X  +  u 

^  (A.l) 

z  =  z  +  w 

and  the  differential  change  in  distance  by  [27] : 


A-1 


dx  =  dx( 1  +  u ,  ) 


dz  =  dx  w , 


(A. 2) 


(ds"^)2 


(dx*)2  +  (dz"*")^ 


By  subtracting  the  first  two  equations  of  Eqn  (A. 2)  into  the 
last  equation  it  follows  that  [27] : 


=  2u.  +  (u.  y  +  (w,  )‘ 


(A. 3) 


Using  the  definition  of  strain  [27] , 


e  =  ( ds  -  dx  )/dx 


(A. 4) 


and  then  rearranging  it  to  get  an  alternate  form  [27] , 


e  +  1/2  e‘ 


(A. 5) 


Substituting  Eqn  (A. 3)  into  Eqn  (A.S)  and  assuming  the 
strain  is  small  so  that  can  be  ignored  gives  [27] : 


e  =  u,  +  l/2(u,"  +  w,‘) 


(A. 6) 


So  far  the  analysis  for  this  line  element  has  been  only 
in  one  plane.  If  this  element  is  in  the  mid-plane  of  a  plate 
element,  the  rotation  of  this  line  element  about  the  normal 
has  to  be  considered.  If  this  rotation  is  included  in  Eqn 
(A. 6)  as  in  Sanders’  nonlinear  shell  equations,  and 
rewritten  to  reflect  mid-plane  strain,  the  result  is: 


e  °  =  u,  +  l/2(u,^  +  w,^  +  <l>^) 

XX  XX 


(A. 7) 


where  in  Sanders’  equations  [3] 


<p  =  l/2(u,  -  V,  ) 

y  X 


(A. 8) 


The  term  that  should  be  included  in  Eqn  (A. 7)  to  represent 
this  rotation  in  the  x-direction  about  the  normal,  for  a 
flat  plate,  is  v,  .  This  does  not  imply  the  use  of  Sander’s 

X 

equation  in  STAGS,  Just  the  importance  of  defining  a  normal 
rotation  term  for  flat  plates  representing  a  shell  geometry. 
Substituting  this  term  in  Eqn  (A. 7)  results  in  [27]: 


=  u,  +  l/2(u,^  +  w,  +  \/.) 


(A. 9) 


Following  this  same  line  of  reasoning,  the  other  two 

mid-plane  terms  e  ®  and  r  °  can  be  written  as  [27] : 

y  *y 


0 

II 

V,  +  l/2(u 
y 

y 

o  ^ 

xy 

u,  +  V,  + 

y  x 

(u,  u 

X 

2 
'  y 


2- 

y ' 


(A. 10) 


Now  if  the  Kirchhoff  hypothesis  is  included,  where  strains 
away  from  the  mid-plane  are  due  to  the  following  curvatures: 


X  =  -w, 

X  XX 


X  =  -W , 

y  yy 


X 

xy 


-w, 

xy 


(A. 11) 


the  full  expression  for  strain  in  the  plate  element  can  be 
written  from  Eqn’s  (A. 9),  (A. 10),  and  (A. 11)  as: 

e  =  e  °  -  zx 

XXX 

e  =  e  ®  -  zx  (A. 12) 

y  y  y 

O 

€  =  €  -  ZX 

xy  xy  xy 

This  kinematic  formulation  allows  for  large  displacements 
and  moderate  rotations  with  small  strains  (due  to  the 
Kirchhoff  hypothesis). 
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Appendix  Bi  Derivation  of  Nonlinear 
Strain  Displacement  Eguat ions 


As  in  STAGS,  one  of  the  primary  strengths  of  SHELL  is 
its  ability  to  analyze  shells  nonlinear ly.  Therefore  a 
short  derivation  of  the  underlying  nonlinear  kinematics  is 
presented.  It  is  important  to  note  that  the  elements  used 
in  SHELL  are  not  plate  elements  but  true  shell  elements. 
Also,  SHELL  incorporates  a  parabolic  distribution  of  through 
the  thickness  shear  strain. 

Therefore,  a  method  to  incorporate  three-dimensional 
(3-D)  effects  is  to  specialize  the  general  3-D  strain 
displacement  relations  expressed  in  arbitrary  orthogonal 
curvilinear  coordinates.  These  are  derived  in  several 
elasticity  texts,  see  Eqn  (2.10)  or  reference  [2]  for 
example. 

The  following  derivation  is  based  on  Linneman’s  [25] 
work.  It  is  presented  here  as  a  means  to  give  the  reader 
the  basic  understanding  of  the  methodology  of  incorporating 
3-D  effects  (transverse  shear)  into  the  2-D  strain 
displacement  relations.  This  derivation  is  not  how  SHELL  is 
formulated,  but  demonstrates  to  the  reader  a  simpler  method 
of  incorporating  transverse  shear  effects  into  the  strain 
displacement  relations  than  is  shown  in  [41] . 

The  coordinate  system  for  a  circular  cylindrical  shell 
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panel  and  the  degrees  of  freedom  used  in  [25]  are  shown  in 
Fig  B.l.  The  x  and  y  axes  are  located  at  the  mid-surface  of 
the  laminate  (z  =  0).  The  degrees  of  freedom  u  (x,y,t), 

O 

V  (x,y,t),  and  w(x,y,t)  are  the  laminate  mid  surface 

O 

displacements  in  the  x,  y,  and  z  directions,  respectively. 
The  degrees  of  freedom  Ui(x,y,t)  and  v>(x,y,t)  are  the 

X  y 

rotations  of  the  laminate  cross  section  from  the  normal  at 
the  mid  surface  with  respect  to  the  x  and  y  axes, 
respectively.  R  is  the  radius  of  curvature,  h  the  laminate 
thickness,  a  the  length  in  the  x-direction  and  b  the  length 


Fig  B.l  Shell  Panel  Coordinates  and  Degrees  of  Freedom 


in  the  y-direction. 

In  order  to  determine  the  displacement  field,  the 
transverse  shear  strains,  y  and  v  ,  need  to  be  modeled. 

»z  yz 

In  classical  laminated  shell  theory,  through  the  thickness 
shear  deformation  is  neglected  according  to  the  Kirchhoff 
-Love  hypothesis  that  plane  cross  sections  remain  plane  and 
perpendicular  to  the  laminate  mid-surface  after  deformation. 
A  displacement  field  that  is  a  first  order  function  of  z  is 
required  in  classical  shell  theory.  Bowlus  [56,57]  and 
Palardy  [58]  in  their  flat  plate  work  modelled  transverse 
shear  strain  using  Mindlin  plate  theory,  which  also  requires 
the  use  of  a  first  order  displacement  field.  Mindlin  plate 
theory  assumes  the  cross  section  remain  plane,  but  is 
allowed  to  rotate  from  the  normal  with  respect  to  the  mid 
surface  after  deformation.  The  assumption  of  no  cross 
sectional  warping  introduces  error,  especially  at  the  top 
and  bottom  surfaces  of  the  laminate,  since  the  model  does 
not  match  the  boundary  conditions  of  zero  transverse  shear 
there.  This  error  is  reduced  by  the  introduction  of  shear 
correction  factor.  SHELL  models  transverse  shear  strain 
parabolically  wherein  the  strains  are  maximum  at  the 
laminate  mid-surface  and  are  zero  at  the  top  and  bottom 
surfaces,  satisfying  the  boundary  conditions.  Fig  B.2 
illustrates  the  various  transverse  shear  concepts  discussed 
above . 
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Plane  Cross  Section 


AFTER  DEFORMATION 


KIrchhoff 


Mindlin 


Parabolic 


Fig  B.2  Transverse  Shear  Strain  Models 
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To  achieve  the  desired  parabolic  transverse  shear,  a 
higher  order  displacement  field  is  required,  as  opposed  to 
the  first  order  displacement  field  used  in  the  Classical  and 
Mindlin  cases.  The  coordinate  displacements  in  the  x  and  y 
directions,  u  and  v.  Mill  be  cubic  functions  of  z:  the 
displacement  in  the  z-direction,  w,  will  be  constant  with 
respect  to  z.  From  Reddy  [59]  and  Saada  [2] ,  the 
generalized  displacement  field  is: 


u(x,y,z,t)  =  u  +  zv  +  z^0  +  z^e 

OX  1  1 

vCx.y.z.t)  = 


(B.l) 


w(x,y,z,t)  =  M 

where  0,0,6,  and  6  will  be  chosen  to  satisfy  the 
12  1  2 

boundary  conditions  of  zero  transverse  shear  strain  at  the 
laminate  top  and  bottom  surfaces. 

Linear  orthogonal  curvilinear  coordinates  are  used  to 
develop  the  strain-displacement  relations  for  simplicity. 
For  a  circular  cylindrical  shell  panel  these  relations 
reduce  to  [25] : 


B-5 


=  U, 


1  + 


xy 


ya 


1  + 


M 

V.  + 

y  R 


u,  +  V, 

y  x 


(B.2) 


1  + 


•"’y  ■  ^ 

V  J 


+  V, 


=  u,  +  w, 

z  X 


where  (  ),  represents  partial  differentiation  with  respect 

to  X  and  so  on.  e  is  assumed  equal  to  zero.  This  implies 

Z 

that  a  change  in  length  in  the  normal  (z)  direction  of  a 

cross  section  perpendicular  to  the  mid  surface  is  not 

considered,  and  is  regarded  as  an  accepted  inconsistency  in 

plate  and  shell  theory.  In  reality,  e  is  not  zero,  but  is 

z 

small  compared  to  the  other  strains.  For  the  laminate,  it 
means  there  are  no  discontinuities  in  e  at  the  lamina 

Z 

boundaries,  but  they  too  are  small. 

The  Donnell  cylindrical  shell  panel  equations  assume 
w  0  in  Eqn  (B.2).  As  shown  in  [25]  ,  this  assumption 


limits  Donnell  theory  to  be  valid  only  for  small 


ratios . 


With  no  transverse  shear,  the  maximum  h/R  limit  under 
Donnell  assumptions  is  approximately  1/500  [59] .  As  will  be 
shown,  with  transverse  shear  included,  the  Donnell  equations 
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are  valid  up  to  h/R  equal  to  approximately  1/50  [59] . 

For  simplicity  assume  at  0  for  the  transverse  shear 
strains,  r  and  r  »  only.  (The  limitations  of  the  model 

xz  yx 

resulting  from  these  simplifications  are  discussed  in 
Appendix  A  of  [25]).  For  the  membrane  strains(  e  ,  e  ,  and 

X  y 

y  )  the  following  polynomial  expansion  is  made: 

xy 


1  + 


Z 

R 


1  - 


z 

R 


(B.3) 


This  approximation  allows  the  strain-displacement  relations 

h 


to  be  valid  for  deep  panels,  with  an 


maximum  limit  of 


approximately  1/5  (See  [24,25]).  The  transverse  shear 
strains  in  Eqn  (B.2)  become: 


7 


y* 


V,  +  w, 

*  y 


R 


7 


XZ 


U,  +  W, 

Z  X 


(B.4) 


If  one  sets  y  (x,y,ih/2,t)  =  O  and  y  (x,y,±h/2,t)  =  0  to 

*yx  ’xz 

satisfy  the  laminate  boundary  conditions,  then  from  Eqn 
(B.l)  it  can  be  shown  that  (see  Appendix  A  [25]): 
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0 


k(v>  +  w,  ) 

X  X 

k(v  +  w,  ) 
y  y 


(B.5) 


where  k  =  - 


The  displacement  field  becomes: 


u(x.y,z,t)  =  u  +  ztj)  +  z  k(v  +  w,  ) 

O  X  X  X 


.y.z.t)  =  [l  *  -|-]v 


+  ZW  +  zk(u>  +  w,  )(B.6) 

o  y  y  y 


w(x,y,z,t)  =  M 


Using  this  displacement  field  in  Eqn  (B.2),  the  strain 
displacement  relations  become: 


u  +  2v  +  z  k(v  +  w,  ) 

o,  X  X,  X  X,  X  XX 


V  +  -5-  +  2V 

o.y  y#y 


+  z^k(v>  +  w,  )  -  -^z^k(v  +  w,  ) 

y.y  yy  R  y.y  yy 


u  +  V  +  z 
o.y  o.x 
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shorthand  notation  can  be  introduced  to  rewrite  the 


strains  as  follows: 
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(Note  the  superscripts  on  the  x  terms  are  not  exponents. 
They  are  for  identification  purposes  only  and  simply 
distinguish  among  the  high  and  low  order  curvature  terms). 
The  strains  at  the  laminate  mid  surface  are  [25] : 
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and  the  curvature  terms  (x)  due  to  bending  and  shear 
deformation  are  defined  as  follows  [25] : 


Appendix  C:  Variational  Formulation 


Both  STAGS  and  SHELL  are  energy  based  finite  element 
programs  dependent  upon  variational  formulation.  The 
following  presentation  of  this  formulation  is  given  for 
better  understanding. 

For  the  static  case  where  the  motion  of  the  system  can 
be  ignored,  the  total  potential  energy  of  a  system,  n  ,  can 

p 

be  given  as: 


n  =  U  -  U 

p 


(C.l  ) 


where  U  is  the  internal  strain  energy  of  the  system  and  U  is 
the  external  work  due  to  applied  forces.  For  a  conservative 
system,  the  change  in  total  potential  energy  is  independent 
of  path  [27] .  The  equations  of  equilibrium  can  then  be 
derived  form  the  first  variation  of  total  potential  energy 
[27].  In  equation  form,  the  principle  of  virtual  work  is: 


60  =  6U  -  5W  =  0 

p 


(C.2) 


The  expression  for  the  internal  strain  energy 


U  = 


f  (e)'^[D](e)dV 

V 


(C.3) 
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where  (e)  is  the  strain  vector  for  the  body  and  [D]  is  the 
material  matrix.  Taking  the  variation  of  Eqn  (C.3)  results 
in: 


6U  = 


£)^[D]6(e)dV 


(C.4) 


or  alternatively  by  taking  the  transpose  of  the  last 
expression  in  Eqn  (C.4)  and  adding  it  to  the  remaining  term 
results  in: 


SU  = 


f  6(e)^[D](£)dV  = 

V 


SC  e)\cf)d\/ 


(C.5) 


Now  one  defines  the  differential  operator ,  [L] ,  which  acts 
on  the  general  displacements,  (u),  to  get  [28]: 

(e)  =  [L](u)  (C.6) 

Next  one  can  define  the  shape  functions,  [N] ,  to  describe 
the  general  displacements  in  terms  of  nodal  displacements 
(parameters),  (a),  of  the  element  where: 


(u)  =  [N](a) 


(C.7) 


If  one  substitutes  Eqn  (C.7)  into  Eqn  (C.6)  the  results 


give: 


[L]  [N](a) 


(C.8) 


(e)  = 

Or  alternatively,  defining  a  matrix  [B]  such  that  [28]: 

[B]  =  [L]  [N]  (C.9) 


Then  Eqn  (C.8)  can  be  written  as: 


(e)  =  [B](a) 


(C.IO) 


Finally,  one  can  take  the  variation  of  Eqn  (C.IO),  transpose 
it,  and  substitute  it  into  Eqn  (C.5)  resulting  in: 


(C.ll  ) 


In  order  to  get  the  full  expression  for  the  variation  of  the 
potential  energy,  assume  that  the  work  term  has  been  defined 
in  terms  of  nodal  displacements,  (a),  and  forces,  (f),  such 
that: 


W  =  (  a  )^(  f  ) 


(C.12) 


Taking  the  variation  of  Eqn  (C.12),  combining  it  with  Eqn 


C-3 


(C.ll),,  and  substituting  these  equations  into  Eqn  (C.2) 
results  in: 


6n  =  6(a)' 

p 


[B]^(a)dV  -  (f)l  =  0 


(C.13) 


Finally,  defining  a  term,  (f),  which  represents  the  sum  of 
the  external  and  internal  forces  (inside  braces,  Eqn  (C.13)) 
[28]  : 


(f)  =  J  [B]’'(<T)dV  -  (f)  = 


or  alternatively  after  integration. 


('F)  =  CK](a)  -  (f)  =  0 


where  [K]  is  the  stiffness  matrix. 


(C.14) 


(C.15) 
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Appendix  D:  Classical  Laminated  Plate  Theory 


The  shells  analyzed  in  this  thesis  are  made  of 
composite  materials.  STAGSC-1 ,  the  computer  program  used  in 
the  analysis,  uses  flat  plates  to  model  a  shell  surface. 
Because  of  this,  a  review  of  composite  plate  theory  is 
presented . 

To  start  the  formulation,  a  few  of  the  classical 
laminated  plate  assumptions  need  to  be  made.  The  laminate 
is  assumed  to  consist  of  perfectly  bonded  plies  where  this 
bond  is  assumed  infinitesimally  thin  and  does  not  allow 
shear  deformation  within  itself  (i.e.  the  displacements 
between  plies  is  continuous:  no  slip  between  plies)  [45] . 
Also,  the  laminate  is  thin  so  that  the  in-plane 
displacements,  u  and  v,  as  well  as  the  normal  displacement, 
w,  are  functions  of  only  two  in-plane  space  variables  (x  and 
y)  [27].  A  plane  stress  assumption  is  assumed  to  where  e  , 

Z 

Y  ,  Y  ,  and  o  are  assumed  to  equal  zero.  Finally,  out  of 

X*  yz  z 

plane  strains  are  due  to  curvatures  as  in  the  Kirchhoff 
hypothesis  for  plates. 

Fig  D.l  shows  the  coordinate  system,  the  associated 
displacements,  and  the  force  and  moment  resultants  (N  ,  M  , 

X  X 

etc.  respectively)  on  the  laminate  [45].  The  x-y  plane  of 
Fig  D.l  coincides  with  the  reference  axis.  From  these 
assumptions,  a  general  expression  can  be  written  for  the 


strains  in  the  laminate  using  the  Kirchhoff  hypothesis  as 
(see  Appendix  A.  Eqn  (A.12)): 
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Mhere  e  ,  e  .and  v  °  represent  in-plane  strains  and  the 

X  y  xy 

terms  x  ,  x  ,  and  x  are  curvatures.  The  curvature  terms 

X  y  xy 

are  defined  as: 
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The  kinematic  relations  for  the  laminate  (Eqn  (D.l))  can  be 
substituted  into  the  constitutive  equation  for  the  laminate 
resulting  in  [45] : 
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where  the  subscript  k  refers  to  the  ply  layer  and  the  CF 
(transformed  reduced  stiffnesses)  quantities  are  defined  as 
follows  [45] : 
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Q  =  Q  cos^e  +  2(Q  +  Q  )sin^0  cos^0  +  Q  sin^0 

11  11  12  66  22 

Q  =(Q  +Q  +4Q  )sin^0  cos^0  +  Q  (sin^0  +  cos^0) 

12  11  22  66  12 

Q  =  Q  sin^0  +  2(Q  +  Q  )sin^0  cos^0  +  Q  cos*d 

22  11  12  66  22 

Q  =  (Q  -  C  -  2Q  )sin  0  cos^0 

16  11  12  66 

+  (Q  -  Q  «  Q  )sin^0  cos  0  (D.4) 

12  22  66 

Q  =  (Q  -  Q  -  2Q  )sin^0  cos  0 

26  11  12  66 

+  (Q  -  Q  +  Q  )sin  0  cos^0 
12  22  66 

Q  =(Q  -Q  -2Q  -4Q  )sin^0  cos^0 

66  11  22  22  66 

+  Q  (sin^0  +  cos^0) 

66 

The  Q  terms  (reduced  stiffnesses)  in  Eqn  (D.4)  are  functions 

of  the  material  properties  as  follows  [45] : 


11 


12 


22 


66 


1  -  y 

12  21 


y  E 
12  2 


U  E 
21  1 


1  -  U  U 

12  21 


1  -  P  P 

12  21 


1  -  P  P 

12  21 


=  G 


12 


(D.5) 


The  angle  0  in  Eqn  (D.4)  is  the  angle  between  the  principle 
axis  (x,y,z)  and  the  material  axis  (1,2,3)  as  in  Fig  D.2. 

The  force  and  moment  resultants  ( forces  and  moments  per 
unit  length)  acting  on  a  composite  plate  can  be  found  by 


D-4 


integrating  the  stresses  in  each  ply  through  the  plate 
thickness,  for  example  [45]: 


N 


M 


pt/2 

J 

-  t/2 

pt/2 

a  z  dz 

-  t/2 


These  force  and  moment  resultants  are  shown  in  Fig  D 
of  the  force  and  moment  resultants  are  put  in  vector 
for  an  iM-layered  laminate  and  the  results  are  [45]  : 
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.1.  All 
form 


(D.7) 
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(a)  Fiber  (1  -2)  Axis  System  (b)  Transformation  into 
Structural  Axis  System 

Figure  D.2  Principal  (x.y.z)  and  Material  (1,2,3) 
Axis  Systems 
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Substituting  Egn  (D.3)  into  Eqns  (D.7)  and  (D.8) 


yields: 
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The  stiffness  matrix  [0]  has  been  removed  from  the  integrals 
in  Eqn  (D.9)  because  it  is  constant  within  a  lamina  but  must 
included  within  the  summation  of  force  and  moment  resultants 
for  each  ply  [45] .  Since  strain  and  curvature  terms  in  Eqn 
(D.9)  are  not  functions  of  z  but  are  mid-plane  values,  they 
can  be  removed  from  the  integration  and  summation  signs 
[45] .  After  integrating  the  remaining  terms,  the  results 
are  [45]  : 
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Appendix  E=  Anisotropic  Cylindrical  Shell 
Theory  Including  Through  the  Thickness 
Shear  Strain 


Lamination  theory  incorporates  constitutive 
relationships  for  orthotropic  lamina  through  the  shell  panel 
thickness  resulting  in  expressions  which  approximate  force 
resultants  in  terms  of  displacement  functions.  This  theory 
provides  concepts  which  are  required  in  the  subsequent 
development  of  the  equations  of  motion  and  boundary 
conditions.  The  constitutive  relationships  are  developed 
for  the  basic  lamina  and  then  to  the  end  result,  the 
structural  laminate.  The  end  results  of  this  section  are 
the  laminate  stiffness  terms  and  force  resultants.  The 
following  derivation  is  based  upon  Appendix  B  and  [25]  for 
simplicity . 

The  plane  stress  constitutive  relationships  for  a 
single  orthotropic  layer  in  the  principle  coordinate  system 
shown  in  Fig  E.l  are  [45]: 
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Note  that  a  =  0  for  plane  stress.  That  is,  the  individual 

Z 

lamina  are  considered  thin  enough  that  the  average  of  a 


Figure  E.l 


Lamina  Material  Coordinates 
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across  the  thickness  is  negligible.  The  S  terms  are 

i  j 

compliance  terms  and  may  be  written  in  terms  of  the  lamina 
engineering  constants  as  [45] : 


(E.2) 


where  E.  are  Young’s  moduli  in  the  i  direction,  u  is 
1  th  ij 

“oisson’s  ration  for  transverse  strain  in  the  j  direction 

th 

when  stressed  in  the  i  direction,  and  G  is  the  shear 

th  i  J 

modulus  in  the  i-j  plane. 

Eqn  (E.l)  may  be  inverted  to  give  the  relationships  of 
the  stresses  in  terms  of  the  strains  [45] : 
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where  Q  are  the  reduced  stiffness  terms  and  are  defined 
i  j 

as  [45]  : 
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A  structural  laminate  consists  of  N  lamina  oriented  at 
different  angles  with  respect  to  each  other .  The  previous 
constitutive  relations  apply  only  to  Fig  E.l  where  the 
lamina-fixed  1-2  axis  system  is  aligned  with  the  laminate 
(or  global)  x-y  axis  system.  If  the  1-2  axis  system  is  not 
aligned  with  the  x-y  axis  system  but  rather  is  at  an  angle 


0  (see  Fig  E.2),  the  reduced  stiffness  matrix,  [Q .  ],  must 
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be  transformed.  The  transformation  matrices  applied  to  the 
stiffness  terms  in  Eqn  (E.3)  to  reflect  the  shift  in  the 
lamina  axes  are  defined  below  [45] ^ 
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where  c  =  cos(0)  and  s  =  sin(0).  The  transformed  stiffness 
matrices  become  [45] : 


Fig  E.2  Arbitrary  Lamina  Coordinates 
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From  [37] ,  the  lamina  constitutive  relationships  can  now  be 
expressed  in  laminate  coordinates  as: 
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where  k  denotes  the  k  lamina  and  the  individual  0  are 

th  i  J 

computed  as  [45] : 
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Finally,  as  an  example,  we  can  substitute  in  the 
expressions  for  the  strains  in  Appendix  B,  Eqn  (B.8)  into 
the  the  constitutive  relations  in  Eqn  (E.7).  The  stress  in 
the  k  lamina  of  the  structural  laminate  is  expressed  as 

t  h 

[45]  : 
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The  resultant  forces  and  moments  and  the  higher  order 
resultant  quantities  acting  on  the  laminate  are  obtained  by 
integrating  the  stresses  in  each  lamina  through  the  laminate 
thickness.  Thus,  for  a  laminate  with  N  lamina  shown  in  Fig 
E.3.  The  resultant  forces  and  moments  and  higher  order 
quantities  are: 
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y*kl.z^)dz 


where  {N  >  and  {Q  }  are  the  resultant  forces,  {M  }  are  the 
i  i  * 

resultant  moments,  and  {S^},  {P^},  {L^},  and  (Rj)  the 

higher  order  quantities  resulting  from  the  higher  order 
strain  displacement  relations. 


Fig  E.3  Geometry  of  an  N-Layered  Laminate  [45] 
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By  substituting  Eqn  (E.9)  into  Eqn  (E.IO),  thereby 
expressing  the  stresses  in  terms  of  the  mid-surface 
displacement  quantities  and  the  transformed  reduced 
stiffness  matrices,  the  integration  is  simplified  because 
the  mid-surface  values  are  independent  of  z  and  can  come  out 
of  the  integral  and  summation  signs  [45] .  This  allows  the 
following  notation  to  be  adopted  for  the  integrated  laminate 
stiffness  matrices: 


i  .  0 

f 

E  , 

F 

i  J 

i  J 

1  J 

1  J 

r  Q 

Q 

Q  1 

11 

12 

16 

Q 

0 

Q 

12 

22 

26 

Q 

0 

Q 

' 

•-  16 

26 

66  -I 

k 

G  ,  H 

ij  i  j 


i  J 


,  .  234567 

(l,z,z  ,z  ,z  ,z  ,z  ,z 


k-1 


1  ,  J 


(E.lla) 

,z®)  dz 

=  1,2,6 


For  the  transverse  shear: 


(A 


IV 

,  D  ,  F  )  =  ) 

iJ  ij  iJ 


k  >  1 


r  Q 

Q  1 

A 

44 

45 

Q 

Q 

45 

55  J 

k 

k-1 


(E.llb) 

i.j  =  4,5 


( 1  ,z^,z^  )dz 


Now,  Eqn  (E.IO)  may  be  written  as: 
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r  > 

N 

1 

N 

2 

► 

K1 

2 

► 

►  = 

KJ 

M 

► 

KJ 

f  'x 

L 

1 

L 

2 

► 

rQ_ 


r«3i 


CA,,]  [B,:  [0  1  [E,.]  [F.,] 

Ij  »J  »J  ij  IJ 


[B..]  [D,.]  [E..]  [F..]  [G.,] 

ij  ij  »J  »J 


[D,]  [E  l  [F.l  [6  ]  [H.,] 

ij  ij  ij  ij  »j 


[E.  .]  [F.  .]  [G.  .]  [H.  ]  [I  ] 

ij  ij  ij  ij  ij 


[F  ]  [G  ]  [H  ]  [I  ]  [J  ] 
ij  ij  ij  ij  ij 


(E.12a) 


xy  J 

i.j  =  1,2,6 
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45 
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55 
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44 


45 


45 


55 


45 


55 


y* 


yz. 


(E.12b) 


where  the  large  matrix  above  is  ( 15  x  15)  and  each  of  its 
sub-matrices  are  the  (3x3)  matrices  in  Eqns  (E.lla)  and 
(E.llb). 
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Appendix  F:  Derivation  of  the  Tangent 
Stiffness  Matrix 


The  tangent  stiffness  matrix  is  a  nonlinear  stiffness 
matrix  used  in  the  Newton-Raphson  (see  section  2.2.7)  for 
solving  a  nonlinear  set  of  equations.  There  are  other 
techniques  for  solving  these  equations,  but  STAGSC-1  and 
SHELL  uses  the  Newton-Raphson  (or  modified  Newton-Raphson), 
and  thus  the  derivation  of  the  tangent  stiffness  matrix  is 
of  interest  in  this  thesis.  This  appendix  derives  the 
tangent  stiffness  matrix  for  STAGS  only.  The  tangent 
stiffness  matrix  development  for  SHELL  is  similar  but  has 
more  terms  and  is  much  more  complex .  See  reference  [37]  for 
a  more  complete  derivation  of  SHELL’S  stiffness  matrix. 

The  tangent  stiffness  matrix  is  derived  for  a  flat 
plate  element,  a  STAGS  type  element,  in  a  general  way 
without  any  reference  to  any  specific  shape  functions.  See 
reference  [43]  for  the  specific  shape  functions  used  in 
STAGS.  Although  the  derivation  is  conducted  in  general 
coordinates,  STAGS  uses  isoparametric  coordinates  in  its 
nonlinear  stiffness  matrix  formulation. 

To  start  the  derivation,  a  form  of  the  tangent 
stiffness  matrix  is  found  from  the  energy  expression. 
Consider  (f),  the  sum  of  external  and  internal  generalized 
forces  (Appendix  C,  Eqn  (C.14)),  which  is  given  as  [28]: 


F-1 


('?)=[  CB]^(a)dV  -(f) 

V 


(F.l) 


where , 

(a)  =  vector  of  stresses 
(f)  =  nodal  forces  of  the  element 
and  [B]  is  defined  as  the  product  of  the  differential 
operator  matrix,  [L] ,  operating  on  the  element  shape 
functir'is,  [N]  .  In  an  equilibrium  state  (S')  (Eqn  (F.l)) 
will  equal  zero.  The  strain  displacement  relations  can  be 
written  (to  include  the  [B]  matrix)  as: 

(e)  =  [L]  [N](a)  =  [B](a)  (F.2) 

where  (a)  is  the  nodal  displacement  vector.  In  the  case  of 
a  nonlinear  stiffness  matrix  (i.e.  tangent  stiffness 
matrix),  [B]  is  redefined  as  [28j : 

[B]  =  [Bo]  +  [Bl]  (F.3) 

where , 

[Bo]  =  constant 

[Bl]  -  function  of  displacements 
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In  order  to  use  the  Newton-Raphson  method,  a  relation 
between  5(a)  and  5('i')  must  be  found.  Taking  the  variation 
of  Eqn  (F.l)  with  respect  to  5(a)  gives  the  relation  needed 
and  is  given  as  [28] : 


5('F) 


[Kt]  5(  a  ) 


(F.4) 


where  [Kx]  is  the  tangent  stiffness  matrix.  In  this 
derivation,  strain  is  assumed  small  therefore  the  equation. 


(a)  = 


[D]  ( e)  -  (  Eo)  +  (  ao) 


(F.5) 


still  applies,  where  [D]  is  the  material  matrix  and  the 
subscripts  indicate  initial  values  of  stresses  and  strains. 
With  Eqns  (F.2),  (F.3),  and  (F.5),  the  variational  terms  in 
Eqn  (F.4)  can  be  rewritten  as: 


5(a)  =  [D]5(£)  =  [D][B]5(a) 


(F.6a) 


and 


5[B]^ 


6[Bl]^ 


(F.6b) 
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With  Eqns  (F.6a)  and  (F.6b),  Eqn  (F.4)  becomes: 


SC'V)  =  J  6[BL]^(a)dV  +  J  [B]  ^  [D]  [B]  dV6(  a  ) 


(F.7) 


The  first  term  in  Eqn  (F.7)  contains  the  initial  stress 
matrix,  [Ka] ,  while  the  last  term  turns  out  to  be  the  linear 
and  the  nonlinear  stiffness  matrices,  [Ko]  and  [Kl] 
respectively,  after  substituting  Eqn  (F.3)  in  for  [B]  [28]. 

The  three  stiffness  matrices  contained  in  Eqn  (F.7)  are: 


[Ko]  = 

r  [Bo] ’’[D]  [Bo]dV 

V 

[Kl]  = 

J  [[Bo]^[D]  [Bl]  +  [Bl]^[D][Bl] 

+  [Bl]^[D]  [BL]jdV 

[Ka]  5(  a  ) 

=  f  5[BL]^(a)dV 

u 

V 


(F.8) 


The  full  expression  for  the  tangent  stiffness  matrix  can  now 
be  written  as: 


[Kt]  =  [Ko]  +  [Kl]  +  [Ka] 


(F.9) 


With  an  expression  for  the  elements  contained  in  the  tangent 
stiffness  matrix  (Eqn  (F.9)),  the  derivation  can  now 


determine  the  terms  contained  in  the  matrices  that  make  up 
the  tangent  stiffness  matrix. 


The  next  steps  in  the  derivation  include  formulating 
the  kinematic  relations,  introducing  the  material  matrix, 
and  giving  the  displacement  relations.  Recalling  the  strain 
expression  from  Appendix  A  ( Eqn  (A.12))  written  in  vector 
form  and  rearranged  to  reflect  in-plane  and  bending  strains 
[28]  : 


(e)  = 


O 

xy 

-w. 

XX 

-w. 

yy 

-2w, 

xy 

(F.IO) 


where  the  distance  for  the  mid-plane,  z,  is  incorporated  in 
the  material  matrix  (see  Appendix  D).  If  Eqn  (F.4)  is 
rewritten,  separating  linear  and  nonlinear  terms  due  to 
in-plane  and  bending  strains,  the  results  are  [28]: 


(e) 


r  A 

£P 

o 

€** 

o 

V  J 


r  A 

oP 


S  i 


(F.ll) 


F-5 


where 


(e"^)  =  linear  in-plane  strains 

O 

(£**)  =  linear  bending  strains 


(e^)  =  nonlinear  in-plane  strains  in  u  and  v 
(e^)  =  nonlinear  in-plane  strains  in  w 
Incorporating  the  linear  constitutive  relationship, 
material  matrix  [D] ,  for  a  composite  material  (derived  in 
Appendix  D),  is  given  as: 


the 


[D] 


[A]  [B] 

[B]  [D] 


(F.12) 


Since  the  distance  from  the  mid-plane,  z,  has  been 
incorporated  in  this  matrix,  the  integration  indicated  in 
Eqn  (F.9)  for  [Ko]  ,  [Kl]  ,  and  [Ka]  reduces  to  an  area 
integral . 

Finally,  the  displacements  are  defined  in  terms  of 
nodal  displacements  using  the  shape  functions  for  the  plate 
element.  If  for  example,  an  element  similar  to  the  QUAF  410 
is  considered  (see  [37] ),  which  has  four  corner  nodes  and 
six  degrees  of  freedom  per  node,  the  displacements  can  be 
given  as: 


CN](a  ) 


(F.13) 
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The  vector  of  element  nodal  displacements,  (a  ),  can  be 

e 

subdivided  into  displacements  that  influence  in-plane 
(superscript  p)  and  bending  (superscript  b)  as  [28]: 


(F.14) 


where  represents  an  average  rotation  about  the  normal 
(similar  to  4>  in  Eqn  (A.8)). 

In  order  to  proceed  further  it  is  necessary  to  expand 
the  expression  for  [B] .  From  Eqn  (F.3)  it  was  shown  that 
[28]  : 


[B]  =  [Bo]  +  [Bl] 


(F.15) 


where , 


[Bo] 


[B**]  O 

(3x?2) 

O 

(3x^2) 


(F  .16a  ) 


and 
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The  matrices  [B**]  (linear  in-plane)  and  [B^]  (linear 

O  O 

bending)  are  standard  matrices  and  are  derived  in  reference 
[28]  with  the  exception  that  in  the  derivation  [B^  the 

O 

normal  rotation  term,  p  (Eqn  (F.14)),  must  be  included  as  a 
degree  of  freedom  at  each  node.  THe  matrices  [Bp  and  [Bp 
are  found  by  taking  a  variation  of  (e^)  and  (e|^)  (E^n 
(F.ll))  respectively,  with  respect  to  the  nodal  degrees  of 
freedom  ( a ) . 

In  matrix  form  (e^)  can  be  written  as: 


(F.17) 


or  [28] : 


(€p  =  I  iRn  (e”) 


(F.18) 


The  vector  (0**)  can  be  related  to  the  in-plane  (u  and  v) 


nodal  degrees  of  freedom  as: 


(e”)  = 


(F.19) 


where  [G'^  is  a  matrix  of  coordinates  and  (a**)  is  the  vector 
of  nodal  displacements  (Eqn  (F.14)).  The  next  step  involves 
taking  the  variation  (e[)-  To  do  this  Eqn  (F.17)  is 
rewritten  as: 


1/2  u,^  +  1/2  v,^ 

1/2  u,^  +  1/2  v,^ 
y  y 


u,^u, 

X  y 


+  V,  V, 

X  y 


(F.20) 


Taking  the  variation  of  Eqn  (F.20)  gives: 


5(c[)  = 


u ,  5u ,  +  V ,  6v , 

XX  XX 

u ,  6u ,  +  V,  5v, 

y  y  y  y 

u,  6u,  +  u,  6u,  +  V,  5v,  +  V,  6v, 

X  y  yx  xy  yx 


(F.21) 


Rewriting  yields: 


«ep  - 


O  u , 


0  u, 

o'  V. 

Vi* 


(F.22) 


where  the  first  matrix  on  the  right  hand  side  equals  [R'*] 


(Eqn  (F.18))  and  the  last  matrix  can  be  written  from  Eqn 
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(F.19)  as: 


u, 

V, 


yj 


(F.23) 


since  [G**]  is  a  function  of  coordinates  only. 

Using  the  definition  of  [R*^  in  Eqn  (F.18)  along  with 
Eqn  (F.23),  Eqn  (F.22)  can  be  written  as: 


6(e[)  =  [r'I  [GP]6(aP)  (F.24) 

From  Eqn  (F.24)  the  resulting  expression  for  [Bp  is,  by 
definition , 


[Bp 

(3x12) 


(F.25) 


Following  the  same  line  of  reasoning  used  to  find  [Bp 
from  (epf  one  can  find  from: 


(F.26) 


or  from  [28] : 
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(£j)  =  1/2  [R^Ke**) 


(F.27) 


that: 


[bJ]  =  [R**]  [G**]  (F.28) 

(3x12)  (3x2)  (2x12) 


where  [G**]  is  composed  of  derivatives  of  the  shape  functions 
that  are  contained  in  the  expression  for  w. 

With  the  expressions  for  [Bo]  and  [Bl]  determined  (Eqns 
(F.16a)  and  (F.16b)),  the  material  matrix  defined  ( Eqn 
(F.12)),  and  recalling  that  the  volume  integral  has  been 
reduced  to  an  area  integral  (Eqn  (F.12)),  the  expressions 
for  linear  and  nonlinear  stiffness  matrices  ( [Ko]  and  [Kl]  ) 
can  be  determined.  If  one  uses  Eqn  (F.8)  and,  after  some 
manipulation , 


.Ksl.  =  J. 


[B'l^CA]  [8*1  [B*1^[B]  [8*1 

o  o  o  o 

[Bl’'[B][B‘1  [B'1^[D][B'1 


dA 


(F.29) 


and 
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• 

p  Cl]  [2] 

[2]  [3] 

(F.30) 

• 

where , 

[1]  = 

[Bf]  t 

o  L 

[B[]^CA][B[]  + 

[Bfl^CA]  [bP] 

L  o 

(F.31 ) 

• 

II 

[Bf]  + 

o  L 

[B[]^[A][B^  + 

[Bfl^CB]  [B**] 

L  o 

(F.32) 

[3]  = 

[B^I^CBlCBh  + 

o  L 

[B^^[A][B^  + 

[Bh^CB]  [B*^ 

L  o 

(F.33) 

• 

The 

final  expression 

necessary  for 

the  tangent 

stiffness  matrix  is  the  initial  stress  matrix  [Ka]  . 
Recalling  Eqn  (F.8),  rewritten  for  convenience  as: 


[Ka]5(a)  =  J  6[Bl]^  (o' )  dV  (F.34) 


The  stresses,  (a'),  are  defined  in  terms  of  the  in-plane 
(superscript  p)  and  bending  (superscript  b)  stresses  as 
[28]  : 


II 

In  N  N  M  M  M  1 

II 

1 

1  X  y  xy  X  y  xyl 

(F.35) 


where  the  primes  on  the  stresses  indicates  stress 
resultants,  not  true  stresses,  and  the  true  stress 


F-12 


resultants  are  given  in  Appendix  D  (Eqn  (D.6)).  The  stress 
resultants  are  average  values  over  the  element.  This  also 
implies  that  integration  through  the  thickness  has  been 
completed  thereby  reducing  the  integration  of  Eqn  (F.34)  to 
that  of  an  area.  Now,  define  the  variation  of  [Bl]  in  Eqn 
(F.34)  from  Eqn  (F.16b)  as: 


■6LBp  O' 

.6[B^  0, 


(F.36) 


Substituting  Eqns  (F.24)  and  (F.25)  into  Eqn  (F.36)  and  then 
substituting  that  result,  as  well  as  Eqn  (F.35),  into  Eqn 
(F.34)  gives  [28]: 

X 

y 

xy 

X 

y 

xy 

Expanding  Eqn  (F.37),  rewriting  it  in  terms  of  in-plane 
(superscript  p)  and  bending  (superscript  b)  expressions  and 
then  expressions  results  in: 
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X 

N 

y 
N 

V  xyj 


>  dA 


(F.38a) 


and 


-  I  MX 


"n  ' 

X 

N 

y 
N 

V.  xyj 


dA 


(F.38b) 


The  following  steps  involve  finding  expressions  for  [Ka**] 
and  [Ka**]  . 

Starting  with  Eqn  (F.38a)  and  rewriting  the  last  two 
matrices  on  the  right  hand  side  using  Eqn  (F.17)  gives: 


• 

r  ■y 

N 

X 

'Ski, 

X 

0 

0 

5u, 

y 

Ski, 

N 

y 

*  M 

fiv. 

y 

0 

X 

5v/» 

N 

^  xyJ 

X 

0 

5v, 

y 

5n/, 

(F.39) 


Expanding  the  right  side  of  Eqn  (F.39)  and  also  taking  the 
variation,  recalling  that  N  ,  N  ,  and  N  are  constants  (Eqn 

X  y  xy 

(F.35)),  results  in: 


F-14 


6[rP]  NN 


N  6u, 

+ 

N 

5u . 

X 

X 

xy 

N  6v, 

+ 

N 

5u. 

y 

y 

xy 

N  6v, 

+ 

N 

6v , 

X 

X 

xy 

N  6v, 

+ 

N 

6v , 

^  y 

y 

xy 

(F.40) 


Rewriting  Eqn  (F.40)  in  matrix  form  gives: 


6[RP]^(N 


V  yJ 


(F.41) 


Recalling  Eqn  (F.23)  and  substituting  into  Eqn  (F.41)  gives: 


[G'’]  ^(aP) 


(F.42) 


Substituting  Eqn  (F.42)  into  Eqn  (F.38a)  results  in  an 
expression  for  [Ka'^  : 


(F.43) 
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With  the  same  type  of  derivation  used  to  find  [Ka**]  , 

[Ka*^  can  be  found.  The  resulting  expression  for  [Ka**]  is: 


''a 


(F.44) 


With  Eqns  (F,43)  and  (F.44),  the  full  expression  for 
the  initial  stress  matrix  can  be  written  as: 


(F.45) 
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Appendix  G: Using  STAGSC-1 


The  STAGSC-1  finite  element  code  was  originally 
developed  for  use  on  the  Cyber.  The  1986  STAGSC-1  version 
is  capable  of  running  both  on  the  Cyber  and  VMS  based 
systems  such  as  Digital  Equipment’s  VAX  11/785  or  VAX  8550. 
Uhen  running  STAGS  on  these  machines,  because  of  the  limited 
time  allowed  for  interactive  use,  each  model  must  be 
submitted  in  batch  form.  Fig  G.l  demonstrates  a  typical 
.COM  file  that  will  run  STAGS  in  batch  form.  Fig  G.2  shows 
the  proper  VMS  syntax  needed  to  submit  the  .COM  file  to  the 
STAGS  queue. 


SUBFILE.COM: 

$  SET  VERIFY 

S  SET  DEFAULT  GAE89D : [SSCHIMME .STAGS] 
S  ASSIGN  FILENAME. INP  F0R005 
S  ASSIGN  FILENAME. OUTl  F0R006 
S  ASSIGN  FILENAME. MOD  F0R021 
$  ASSIGN  FILENAME.DAT  F0R002 

•  ASSIGN  FILENAME. INP  F0R025 
$  RUN  STAGS$DIR:STAGS1 

t  COPY  FILENAME.DAT  FILENAME.DAT 
S  ASSIGN  FILENAME.BIN  F0R005 

•  ASSIGN  FILENAME. 0UT2  F0R006 
S  ASSIGN  FILENAME. RST  F0R020 
S  ASSIGN  FILENAME. SOD  F0R022 
S  RUN  STAGS«DIR:ET 

S  DEL  F0R009.DAT 
S  PUR  SUBFILE.COM 
«  EXIT 


Fig  G.l 


Typical  Command  File  for  Batch  Execution 


$  SUBMIT/QUE=SYS^STAGS/N0TIFY/NAME=STA6S/KEEP/L0G=GAE89D : 
[SSCHIMME .STAGS] FILENAME .LOG/NOPRINTER  SUBFILE .COM 


Fig  G.2  VMS  Syntax  for  Submitting  Batch  File 

Unlike  the  STAGS  Cyber  version,  the  VMS  based  code 
requires  two  separate  files  for  STAGS  to  run.  The  first 
file  is  called  the  .INP  file.  It  is  exactly  the  same  as  the 
input  deck  file  for  the  Cyber  based  code.  The  second  file 
is  called  the  .BIN  file.  This  file  controls  the  STAGS2 
phase  for  the  solution  algorithm.  Examples  of  both  files 
are  given  in  Figs  G.3  through  G.IO.  The  reader  is  suggested 
to  refer  heavily  to  Ref  27  while  looking  at  these  examples 
for  clarification. 

Some  of  the  models  used  for  STAGS  required  tremendously 
large  amount  of  CPU  time  ( 100000+  CPU  seconds ) .  On  the  D1 
record,  the  user  specifies  the  amount  of  CPU  time  STAGS  is 
supposed  to  run  for  the  model.  The  author  normally  used 
90000  (because  of  the  combination  of  using  0.5*  mesh  and 
using  the  QUAF  411  element)  as  most  of  the  models  considered 
large  rotations  (m  30**).  These  large  rotations  required  a 
substantial  number  iterations  and  cutting  the  load  step  in 
order  to  converge  to  a  solution.  This  tends  to  drive  the 
required  CPU  time  for  solution  upward  in  a  somewhat 
exponential  manner.  STAGS  allows  for  only  5  characters  for 


S121.INP; 


NONLINEAR  ANALYSIS,  12“  PANEL.  NO  CUTOUT  $A1 

3.1  ,1 ,1 ,0,0,1  $B1 

1  $B2 

1,0,1  $B3 

.0001 , .0002, .04  $C1 

•  0,90000,40,-20,-1 , .00001  $01 

25,25  $F1 

1  $11 

20. 461E6, 0.0205, .8638E6,1 . ,1 . ,1.3404E6,1 .  $12 

1,1,8  $K1 

1..  005. 0.0.1  $K2 

•  1 , .005,-45.0,0  $K2 

1 . . 005.45.0.0  $K2 

1 . . 005.90.0.1  $K2 

1 . . 005.90.0.1  $K2 

1 . . 005.45.0.1  $K2 

1 . . 005,-45.0,1  $K2 

•  1.. 005, 0.0.1  $K2 

5  $M1 

0.0,12.0,0.0,57.2958,12.0  $M2A 

1  $M5 

411  $N1 

0,0, 0,0,0  $P1 

•  100,000  $P2 

100,100  $P2 

000,000  $P2 

100,100  $P2 

1  $Q1 

1.1  $02 

•  1.0, -1,1, 1,0  $03 

1,0, 0,0, 0,1  $R1 


S121.BIN 


NONLINEAR  ANALYSIS,  12“  PANEL,  NO  CUTOUT 
.0001 , .0005, .04 
0 ,90000,40,-20 ,-l , .00001 


$P1 

$P2 

$P3 

$P4 


Fig  G.3  Input  Deck  for  12x12  Panel,  Simply  Supported, 
No  Cutout 


G-3 


S122.INP: 


NONLINEAR  ANALYSIS,  12"  PANEL,  4“  CUTOUT,  CENTERED  SAl 

3.1 .1 .1  ,0,0,1  $B1 

1  $B2 

1,0,1  $B3 

.0001 , .0002, .04  SCI 

0,90000,40,-20,-1 , .00001  SDl 

25,25  SFl 

1  SIl 

20.461E6,0.0205, .8638E6,1. .1. ,1.3404E6,1.  SI2 

1,1,8  SKI 

1. . 005.0.0.1  SK2 

1 . . 005,-45.0,0  SK2 

1 . . 005.45.0.0  SK2 

1,-005,90.0,1  SK2 

1 . . 005.90.0.1  SK2 

1 . . 005.45.0.1  SK2 

1 . . 005,-45.0,1  SK2 

1..  005. 0.0.1  SK2 

5  SMI 

0.0,12.0,0.0,57.2958,12.0  ‘  SM2A 

1  SM5 

411,0,0,1  SNl 

9,17,9,17  SN8 

0,0, 0,0,0  SPl 

100,000  SP2 

100,100  SP2 

000,000  SP2 

100,100  SP2 

1  SOI 

1 . 1  SQ2 

1.0, -1,1, 1,0  SQ3 

1,0, 0,0, 0,1  SRI 


S122.BIN 

NONLINEAR  ANALYSIS,  12"  PANEL,  4"  CUTOUT,  CENTERED  SPl 
.0001 , .0005, .04  SP2 
0,90000,40,-20,-1,-00001  SP3 
1  SP4 


Fig  G.4  Input  Deck  for  12x12  Panel,  Simply  Supported, 
4"  Cutout,  Centered 


S123.INP: 


NONLINEAR  ANALYSIS,  12 •*  PANEL,  4"  CUTOUT,  1“  OFFSET  SAl 
3 , 1 , 1 , 1 ,0 ,0 , 1  $B1 

1  $B2 

1,0,1  $B3 

.0001 , .0002, .04  $C1 

0 ,90000 ,40 ,-20 ,-l , .00001  $D1 

25,25  $F1 

1  SIl 

20. 461E6, 0.0205, .8638E6,1. ,1. ,1.3404E6,1.  SI2 

1.1,8  $K1 

1..  005. 0.0.1  $K2 

1 . . 005 ,-45 .0  ,0  $K2 

1 . . 005 .45 .0 .0  $K2 

1. . 005.90.0.1  $K2 

1. . 005.90.0.1  SK2 

1 . . 005.45.0.1  $K2 

1 . . 005 ,-45 .0 ,1  SK2 

1..  005. 0.0.1  $K2 

5  $M1 

0.0,12.0,0.0,57.2958,12.0  SM2A 

1  ^M5 

411,0,0,1  $N1 

9,17,7,15  $N8 

0,0,0 ,0,0  ^P1 

100,000  $P2 

100,100  $P2 

000,000  *P2 

100,100  SP2 

1  $Q1 

1.1  $Q2 

1.0,— 1,1, 1,0  $Q3 

1,0, 0,0, 0,1  SRI 


S121.BIN 

NONLINEAR  ANALYSIS,  12"  PANEL,  4"  CUTOUT,  1"  OFFSET  SPl 
.0001 , .0005, .04  SP2 
0,90000,40,-20,-1 , .00001  $P3 
1  SP4 


Fig  G.5  Input  Deck  for  12x12  Panel,  Simply  Supported, 
4"  Cutout,  1"  Offset 


S82 . INP : 


NONLINEAR  ANALYSIS,  8"  PANEL.  4"  CUTOUT  SAl 

3, 1,1, 1,0, 0,1  $81 

1  $82 

1,0,1  $83 

.0001 , .0002, .04  $C1 

•  0,90000,40,-20,-1 , .00001  $01 

25,17  $F1 

1  $11 

20. 461E6, 0.0205, .8638E6,1 . ,1 . ,1 .3404E6,1 .  $12 

1,1,8  $K1 

1..  005. 0.0.1  $K2 

•  1 , .005,-45.0,0  $K2 

1 . . 005.45.0.0  $K2 

1 . . 005.90.0.1  $K2 

1 . . 005.90.0.1  $K2 

1 . . 005.45.0.1  $K2 

1 . . 005,-45.0,1  $K2 

•  1,. 005, 0.0,1  $K2 

5  $M1 

0.0,12.0,0.0,38.1972,12.0  $M2A 

1  $M5 

411,0,0,1  $N1 

9,17,5,13  $N8 

•  0,0, 0,0,0  $P1 

100,000  $P2 

100,100  $P2 

000,000  $P2 

100,100  $P2 

1  $Q1 

•  1,1  $Q2 

1.0, -1,1, 1,0  $Q3 

1,0, 0,0, 0,1  $R1 


S82.BIN 


NONLINEAR  ANALYSIS,  8'  PANEL,  4"  CUTOUT 
.0001, .0005, .04 
0 ,90000 ,40,-20,-! , .00001 
1 


$P1 

$P2 

$P3 

$P4 


Fig  G.6 


Input  Deck  for  8x12  Panel,  Simply  Supported, 
4"  Cutout,  Centered 


6-6 


F122.INP: 


NONLINEAR  ANALYSIS,  12"  PANEL.  4"  CUTOUT.  FREE  EDGES  $A1 
3, 1,1, 1,0, 0,1  SBl 

1  SB2 

1.0,1  SB3 

.0001 , .0002, .04  SCI 

0,90000,40,-20,-1 , .00001  SDl 

23 , 25  SF 1 

1  $11 

20. 461E6, 0.0205, .8638E6,1 . ,1 . ,1 .3404E6,1  .  SI2 

1,1,8  SKI 

1..  005. 0.0.1  SK2 

1 . . 005 ,-45 .0 ,0  SK2 

1 . . 005 .45 .0 .0  SK2 

1.-005,90. 0,1  SK2 

1 . . 005.90.0.1  SK2 

1. . 005.45.0.1  SK2 

1. . 005,-45.0.1  SK2 

1..  005. 0.0.1  SK2 

5  SMI 

0.0,11.0,0.0,57.2958,12.0  SM2A 

1  $M5 

411,0,0,1  SNl 

8,16,9,17  SN8 

0,0 ,0,0,0  SPl 

100,000  SP2 

111,111  SP2 

000,000  SP2 

111,111  SP2 

1  $Q1 

1 , 1  $02 

1.0, -1,1, 1,0  SQ3 

1,0, 0,0, 0,1  SRI 


F122.BIN 

NONLINEAR  ANALYSIS,  12’  PANEL,  4’  CUTOUT,  CENTERED  $P1 
. 0001 , . 0005 , . 04  SP2 
0,90000,40,-20,-1 , .00001  SP3 
1  $P4 


Fig  G.7  Input  Deck  for  12x11  Panel,  Unsupported, 
4’  Cutout,  Centered 


G-7 


F82.INP 


NONLINEAR  ANALYSIS,  8"  PANEL.  4"  CUTOUT,  CENTERED 
3,1 ,1 ,1 ,0,0,1 
1 

1.0.1 

.0001 , .0002, .04 
0,90000,40,-20,-1 . .00001 
23.17 
1 

20 . 461E6 , 0 . 0205 , . 8638E6 . 1 . , 1 . . 1 . 3404E6 . 1 . 

1.1,8 

1 ,  .005,0.0,1 
1 ,  .005,-45.0,0 
1 ,  .005,45.0,0 
1  ,  .005,90.0,1 

1. . 005.90.0.1 

1 . . 005.45.0.1 

1 . . 005,-45.0,1 

1 . . 005.0.0.1 
5 

0.0,11 .0,0.0,38.1972,12.0 
1 

411 ,0,0,1 
8,16.5,13 
0,0, 0,0,0 
100,000 
111,111 
000,000 
111,111 
1 

1.1 

1.0, -1,1, 1,0 
1.0, 0,0. 0,1 


$A1 

SBl 

$B2 

SB3 

$C1 

$01 

$F1 

$11 

$12 

$K1 

$K2 

$K2 

$K2 

$K2 

$K2 

$K2 

$K2 

$K2 

$M1 

$M2A 

$M5 

$N1 

$N8 

$P1 

$P2 

$P2 

$P2 

$P2 

$01 

$02 

$03 

$R1 


F82.BIN 


NONLINEAR  ANALYSIS,  8“  PANEL,  4* 
.0001 , .0005, .04 
0 ,90000,40 ,-20 ,-l , .00001 
1 


CUTOUT,  CENTERED  $P1 

$P2 

$P3 

$P4 


G.8  Input.  Deck  for  8x11  Panel,  Unsupported, 
4"  Cutout,  Centered 


G-8 


S121-I.INP: 


NONLINEAR  ANALYSIS,  12"  PANEL,  NO  CUTOUT  $A1 

3, 1,1, 1,0, 0,1  $B1 

1  $B2 

1,0,1  SB3 

.0001 , .0002, .04  $C1 

0,90000,40,-20,-1 , .00001  $D1 

25,25  SFl 

1  SIl 

20.461E6,0.0205, .8638E6,1. ,1. ,1.3404E6,1 .  SI2 

1,1,8  $K1 

1..  005. 0.0.1  $K2 

1 . . 005,-45.0,0  SK2 

1 . . 005.45.0.0  SK2 

1. . 005.90.0.1  $K2 

1,-005,90.0,1  $K2 

1 . . 005.45.0.1  $K2 

1 . . 005,-45.0,1  SK2 

1..  005. 0.0.1  SK2 

5  SMI 

0.0,12.0,0.0,57.2958,12.0  SM2A 

1 , 1  SM5 

6.0,28.6479,4.0,19.0986, .005  SM6 

411  SNl 

0,0, 0,0,0  SPl 

100,000  SP2 

100,100  SP2 

000,000  SP2 

100,100  SP2 

1  SQl 

1 , 1  SQ2 

1.0, -1,1, 1,0  SQ3 

1,0, 0,0, 0,1  SRI 


S121I.BIN 


NONLINEAR  ANALYSIS,  12"  PANEL,  NO  CUTOUT  SPl 
. 0001 , . 0005 , . 04  SP2 
0,90000,40,-20,-1, .00001  SP3 
1  SP4 


Fig  G.9  Input  Deck  for  12x12  Panel,  Simply  Supported, 
No  Cutout,  with  Imperfection 


S121V.INP: 


NONLINEAR  ANALYSIS,  12"  PANEL,  NO  CUTOUT  $A1 

3,1 ,1 ,1 ,0,0,1  $B1 

1  SB2 

1,0,1  $B3 

.0001 , .0002, .04  $C1 

0,90000,40,-20,-1 , .00001  $D1 

25,25  $F1 

1  SIl 

20.461E6,0.0205, .8638E6,1 . ,1 . ,1 .3404E6,1  .  SI2 

1,1,8  9K1 

1..  005. 0.0.1  $K2 

1 . . 005,-45.0,0  SK2 

1 . . 005.45.0.0  «K2 

1 . . 005.90.0.1  SK2 

1 . . 005.90.0.1  »K2 

1 . . 005.45.0.1  $K2 

1 . . 005,-45.0,1  SK2 

1..  005. 0.0.1  ♦K2 

5  $M1 

0.0,12,0,0.0,5757.2958,12.0  SM2A 

1  $M5 

411  $N1 

0,0, 0,0,0  $P1 

100,000  $P2 

110,100  $P2 

000,000  $P2 

110,100  SP2 

1  $Q1 

1 , 1  $02 

1.0, -1,1, 1,0  $03 

1,0, 0,0, 0,1  $R1 


S121V.BIN 


NONLINEAR  ANALYSIS,  12"  PANEL,  NO  CUTOUT  $P1 

.0001 , .0005, .04  $P2 
0 ,90000 ,40 ,-20 ,-l , .00001  $P3 
1  $P4 


Fig  G.IO  Input  Deck  for  12x12  Panel,  Simply  Supported, 
No  Cutout,  with  V  Free 


NSEC  variable.  Thus,  the  user  cannot  specify  more  than 
99999  CPU  seconds.  In  this  case,  the  STAGS  restart 
capability  must  be  used.  Fig  G.ll  demonstrates  the  .COM 
file  needed  to  accomplish  a  STAGS  restart. 


SUBFILE.COM: 


S  SET  VERIFY 

•  SET  DEFAULT  GAE89D: [SSCHIMME .STAGS] 
$  ASSIGN  FILENAME.BIN  F0R005 

S  ASSIGN  FILENAME. 0UT3  F0R006 
S  COPY  FILENAME. SOD  FILENAME. RST 
S  COPY  FILENAME.DAT  F0R002 

•  ASSIGN  FILENAME. RST  F0R020 
S  ASSIGN  FILENAME .SOD  F0R022 
$  RUN  STAGS«DIR:ET 

«  PUR  SUBFILE.COM 

•  EXIT 


Fig  G.ll  Typical  Command  File  for  Batch  Restart 

Besides  a  new  .COM  file,  the  .BIN  file  needs  to  be 
modified.  In  the  Cl  deck,  the  starting  load  is  set  to 
correspond  to  the  final  load  completed  in  the  prior  run.  It 
is  strongly  recommended  that  the  user  chooses  not  the  last 
load  step  but  the  second  to  the  last  step.  This  ensures 
that  the  model,  when  restarted,  will  essentially  follow  the 
same  equilibrium  path  before  STAGS  stopped  execution.  The 
increment  should  also  be  matched  to  a  more  appropiate  value. 


usually  the  increment  last  being  used  by  STAGS  before 
execution  was  terminated.  The  step  number  in  the  D1  deck 
changed  to  match  the  incremental  load  value  being  used  as  a 
starting  point  for  the  restart  in  deck  Cl.  Fig  G.12 
demonstrates  a  .BIN  file  modified  for  a  STAGS  restart. 
Submit  the  modified  .BIN  file  through  the  .COM  file  as 
before  (See  Fig  G.2). 


S122.BIN: 


NONLINEAR  ANALYSIS,  12"  PANEL,  4"  CUTOUT,  RESTART  $P1 
3 .9217E-04 . ,00001 , .04  $P2 
27,75000,40,-20,-1, .00001  $P3 
1  *P4 


Fig  G.12 


Modified  Control  File  for  Restart  Execution 


G- 


Appendix  H:  Miscellaneous  Information 


The  folloMing  are  tabulated  results  for  the  STAGS  and 
SHELL  analytical  models  and  the  experimental  results. 

NOTE:  SS  =  Simply  Supported,  US  =  Unsupported,  NC  = 
Cutout ,  CC  =  Centered  Cutout ,  OC  =  Of f center  Cutout . 


Table  H.l  Individual  Panel  Results,  STAGS 
Collapse  Load,  lbs 


Panels 

Perfect 

Imperfect 

V  Free 

12" ,  SS,  NC 

5060 

4950 

4754 

12"  SS,  CC 

2521 

N/A 

N/A 

12"  SS,  OC 

2561 

N/A 

N/A 

8",  SS,  NC 

3761 

3700 

3174 

8",  SS,  CC 

2178 

N/A 

N/A 

8"  ,  SS,  OC 

2465 

N/A 

N/A 

12"  ,  US,  NC 

2298 

2260 

N/A 

12",  US,  CC 

1186 

N/A 

N/A 

12",  US,  OC 

1284 

N/A 

N/A 

8",  US,  NC 

1382 

1350 

N/A 

8",  US,  CC 

323 

N/A 

N/A 

8",  US,  OC 

460 

N/A 

N/A 

Table  H.2  Individual  Panel  Results,  STAGS 


Top  Edge  Displacment ,  in 


Panels 

Perfect 

Imperfect 

V  Free 

12",  SS,  NC 

.016433 

.01842 

.015437 

12"  SS,  CC 

.014823 

N/A 

N/A 

12"  SS,  OC 

.014535 

N/A 

N/A 

8",  SS,  NC 

.018334 

.019454 

.015457 

8",  SS,  CC 

.021156 

N/A 

N/A 

8"  ,  SS,  OC 

.026007 

N/A 

N/A 

12",  US,  NC 

.0068329 

.011109 

N/A 

12",  US,  CC 

.0088393 

N/A 

N/A 

12" ,  US,  OC 

.0087854 

N/A 

N/A 

8",  US,  NC 

.0080207 

.027807 

N/A 

8",  US,  CC 

.0087811 

N/A 

N/A 

8",  US,  OC 


0071766 


N/A 


N/A 


Table  H.3  Individual  Panel  Results,  SHELL 
Collapse  Load,  lbs 


Panels 

Perfect 

12" ,  US,  NC 

12".  US,  CC 

12",  US,  OC 

8 " ,  US ,  NC 

8 "  .  US ,  CC 

312.8 

8",  US,  OC 

Table  H.4 


Individual  Panel  Results, 


Top  Edge  Displacment ,  in 


Panels 

Perfect 

12"  ,  US,  NC 

12" ,  US,  CC 

12" ,  US,  OC 

8",  US,  NC 

8"  ,  US,  CC 

.0098 

8 " ,  US ,  OC 

SHELL 


Table  H.5  Individual  Panel  Results,  Experimental 


Panels 

Collapse  Load,  lb 

Displacment,  in 

12“ ,  SS,  NC 

4851 

.0219 

12“  SS,  CC 

2392 

.0212 

12“  SS,  OC 

2406 

.016 

8  “  ,  SS ,  NC 

3473 

.0242 

8“,  SS,  CC 

1918 

.02165 

8  “  ,  SS ,  OC 

1619 

.022 

12“ ,  US,  NC 

1890 

.0291 

12“ ,  US,  CC 

1093 

.0172 

12“ ,  US,  OC 

1127 

.0182 

8“  ,  US,  NC 

1311 

.0279 

8“,  US,  CC 

300 

.0123 

8 “ ,  US ,  OC 

403 

.0099 

VITA 
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the  vertical  supports  reduced  the  total  collapse  load  by  at  least  507,. 
Surface  imperfections  were  seen  to  reduce  the  collapse  load,  but  more 
importantly  these  imperfections  accurately  predicted  the  nonlinear 
response  of  the  shell.  STAGSC-1  accurately  predicted  the  composite 
panels'  response  to  axial  compression.  Transverse  ^ear  was  seen  to  be 
important  if  the  panels*  free  edges  underwent  15-17^  or  more  of  rotation. 
Hence  the  SHELL  program  gave  better  resul^  than  STAGSC-1  for  those 
panels  that  saw  rotations  greater  than  l/fT  The  residual  stress  test 
showed  there  was  a  measureable  change  in ^strain  due  to  cutting  the  panel. 
However,  the  magnitude  of  this  cha^e  is  minimal  when  compared  to  the 
total  collapse  load  of  the  shell.  < 


